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Material  Chemistry  Structure  C-S-H  Jennit  -  A  Computational  Modeling  Study 

ABSTRACT 

Material  properties  and  performance  are  governed  by  material  molecular  chemistry  structures  and  molecular  level 
interactions.  Methods  to  understand  relationships  between  the  material  properties  and  performance  and  their 
correlation  to  the  molecular  level  chemistry  and  morphology,  and  thus  find  ways  of  manipulating  and  adjusting 
matters  at  the  atomistic  level  in  order  to  improve  material  performance  are  required.  A  computational  material 
modeling  methodology  is  investigated  and  demonstrated  for  a  key  cement  hydrated  component  material  chemistry 
structure  of  Calcium-Silicate-Hydrate  (C-S-H)  Jennite  in  this  work. 

The  effect  of  material  ion  exchanges  on  the  mechanical  stiffness  properties  and  shear  deformation  behavior  of 
hydrated  cement  material  chemistry  structure  of  Calcium  Silicate  Hydrate  (C-S-H)  Jennite  was  studied.  Calcium  ions 
were  replaced  with  Magnesium  ions  in  Jennite  structure  of  the  C-S-H  gel.  Different  level  of  substitution  of  the  ions 
was  used.  The  traditional  Jennite  structure  was  obtained  from  the  American  Mineralogist  Crystal  Structure  Database 
and  super  cells  of  the  structures  were  created  using  a  Molecular  Dynamics  Analyzer  and  Visualizer  Material  Studio. 
Molecular  dynamics  parameters  used  in  the  modeling  analysis  were  determined  by  carrying  out  initial  dynamic 
studies.  64  unit  cell  of  C-S-H  Jennite  was  used  in  material  modeling  analysis  studies  based  on  convergence  results 
obtained  from  the  elastic  modulus  and  total  energies.  NVT  forcite  dynamics  using  compass  force  field  based  on  200 
ps  dynamics  time  was  used  to  determine  mechanical  modulus  of  the  traditional  C-S-H  gel  and  the  Magnesium  ion 
modified  structures.  NVT  Discover  dynamics  and  COMPASS  forcefield  was  used  in  the  material  modeling  studies  to 
investigate  the  influence  of  ionic  exchange  on  the  shear  deformation  of  the  associated  material  chemistry  structures. 
A  prior  established  quasi-static  deformation  method  to  emulate  shear  deformation  of  C-S-H  material  chemistry 
structure  that  is  based  on  a  triclinic  crystal  structure  was  used,  by  deforming  the  triclinic  crystal  structure  at  0.2 
degree  per  time  step  for  75  steps  of  deformation. 

It  was  observed  that  there  is  a  decrease  in  the  total  energies  of  the  systems  as  the  percentage  of  magnesium  ion 
increases  in  the  C-S-H  Jennite  molecular  structure  systems.  Investigation  of  effect  of  ion  exchange  on  the  elastic 
modulus  shows  that  the  elastic  stiffness  modulus  tends  to  decrease  as  the  amount  of  Mg  in  the  systems  increases, 
using  either  compass  or  universal  force  field.  On  the  other  hand,  shear  moduli  obtained  after  deforming  the  structures 
computed  from  the  stress-strain  curve  obtained  from  material  modeling  increases  as  the  amount  of  Mg  increases  in 
the  system.  The  present  investigations  also  showed  that  ultimate  shear  stress  obtained  from  predicted  shear  stress  - 
strain  also  increases  with  amount  of  Mg  in  the  chemistry  structure.  Present  study  clearly  demonstrates  that 
computational  material  modeling  following  molecular  dynamics  analysis  methodology  is  an  effective  way  to  predict 
and  understand  the  effective  material  chemistry,  and  additive  changes  on  the  stiffness  and  deformation  characteristics 
in  cementitious  materials  and  extendable  to  others. 
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Abstract 

Material  properties  and  performance  are  governed  by  material  molecular  chemistry 
structures  and  molecular  level  interactions.  Methods  to  understand  relationships  between  the 
material  properties  and  performance  and  their  correlation  to  the  molecular  level  chemistry  and 
morphology,  and  thus  find  ways  of  manipulating  and  adjusting  matters  at  the  atomistic  level  in 
order  to  improve  material  performance  are  required.  A  computational  material  modeling 
methodology  is  investigated  and  demonstrated  for  a  key  cement  hydrated  component  material 
chemistry  structure  of  Calcium-Silicate-Hydrate  (C-S-H)  Jennite  in  this  work. 

The  effect  of  material  ion  exchanges  on  the  mechanical  stiffness  properties  and  shear 
deformation  behavior  of  hydrated  cement  material  chemistry  structure  of  Calcium  Silicate 
Hydrate  (C-S-H)  Jennite  was  studied.  Calcium  ions  were  replaced  with  Magnesium  ions  in 
Jennite  structure  of  the  C-S-H  gel.  Different  level  of  substitution  of  the  ions  was  used.  The 
traditional  Jennite  structure  was  obtained  from  the  American  Mineralogist  Crystal  Structure 
Database  and  super  cells  of  the  structures  were  created  using  a  Molecular  Dynamics  Analyzer 
and  Visualizer  Material  Studio.  Molecular  dynamics  parameters  used  in  the  modeling  analysis 
were  determined  by  carrying  out  initial  dynamic  studies.  64  unit  cell  of  C-S-H  Jennite  was  used 
in  material  modeling  analysis  studies  based  on  convergence  results  obtained  from  the  elastic 
modulus  and  total  energies.  NVT  forcite  dynamics  using  COMPASS  force  field  based  on  200  ps 
dynamics  time  was  used  to  determine  mechanical  modulus  of  the  traditional  C-S-H  gel  and  the 
Magnesium  ion  modified  structures.  NVT  Discover  dynamics  and  COMPASS  forcefield  was 
used  in  the  material  modeling  studies  to  investigate  the  influence  of  ionic  exchange  on  the  shear 
deformation  of  the  associated  material  chemistry  structures.  A  prior  established  quasi-static 
deformation  method  to  emulate  shear  deformation  of  C-S-H  material  chemistry  structure  that  is 
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based  on  a  triclinic  crystal  structure  was  used,  by  deforming  the  triclinic  crystal  structure  at  0.2 
degree  per  time  step  for  75  steps  of  deformation. 

It  was  observed  that  there  is  a  decrease  in  the  total  energies  of  the  systems  as  the 
percentage  of  magnesium  ion  increases  in  the  C-S-H  Jennite  molecular  structure  systems. 
Investigation  of  effect  of  ion  exchange  on  the  elastic  modulus  shows  that  the  elastic  stiffness 
modulus  tends  to  decrease  as  the  amount  of  Mg  in  the  systems  increases,  using  either  compass  or 
universal  force  field.  On  the  other  hand,  shear  moduli  obtained  after  deforming  the  structures 
computed  from  the  stress-strain  curve  obtained  from  material  modeling  increases  as  the  amount 
of  Mg  increases  in  the  system.  The  present  investigations  also  showed  that  ultimate  shear  stress 
obtained  from  predicted  shear  stress  -  strain  also  increases  with  amount  of  Mg  in  the  chemistry 
structure.  Present  study  clearly  demonstrates  that  computational  material  modeling  following 
molecular  dynamics  analysis  methodology  is  an  effective  way  to  predict  and  understand  the 
effective  material  chemistry,  and  additive  changes  on  the  stiffness  and  deformation 


characteristics  in  cementitious  materials  and  extendable  to  others. 
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CHAPTER  1 
Introduction 

Research  has  suggested  that  at  least  20%  of  all  product  innovation  is  based  on  the 
introduction  of  new  and  innovative  materials.  It  has  also  been  established  that  material  properties 
and  performance  is  governed  by  molecular  structures  and  molecular  level  interactions  [1].  Hence 
it  behooves  engineers  and  scientists  to  find  ways  of  understanding  the  relationships  between  the 
material  properties  and  performance  that  are  influenced  by  molecular  level  chemistry  and 
morphology,  and  thus  find  innovative  ways  of  manipulating  and  adjusting  matters  at  the 
atomistic/molecular  level  in  order  to  improve  material  performances.  For  example,  cementitious 
materials  that  are  the  binders  for  mortar  and  concrete  have  material,  geometrical  features  ranging 
from  material  chemistry/molecular/nano,  micro,  meso,  and  macro  scales.  Features  and  changes  at 
the  material  chemistry/nanoscale  level  influence  the  hydration  process,  formed  micro  scale 
morphology,  associated  properties  and  behavior  at  engineering  length  scales.  In  order  to  engineer 
the  material,  it  is  required  to  develop  a  thorough  understanding  of  the  material  chemistry  and 
processes  that  define  the  material  structure,  this  is  particularly  important  for  complex  material 
systems  such  as  cement.  A  brief  background  discussion  of  the  material  chemistry  of  cement  and 
its  multi-scale  features  is  presented  next. 

1.1  Cement:  Evolution  of  Material  Chemistry 

Because  of  its  chemical  durability,  remarkable  mechanical  properties  and  high  versatility, 
cement  has  remained  the  most  widely  utilized  material  in  the  world  [2],  The  starting  material  of 
cement  is  the  clinker  phase  which  is  the  unhydrated  cement.  The  unhydrated  cement  phase  is 
comprised  of  tri-calcium  silicate  (alite),  di-calcium  silicate  (belite),  tri-calcium  aluminate,  tetra- 
calcium  aluminoferrite,  and  additional  trace  compounds  [3].  Table  1  below  shows  the  various 
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cement  clinker  constituents  and  their  corresponding  mass  percentages  in  a  typical  Portland 
cement  clinker. 

Table  1 


Unhydrated  Cement  Constituent 


Cement  Clinker 

Chemical 

Formula 

Oxide  Formula 

Shorthand 

Notation 

Mass 

% 

Tri-Calcium  Silicate  (alite) 

Ca3Si05 

(Ca0)3Si02 

C3S 

50-70 

Di-Calcium  Silicate  (belite) 

Ca2Si04 

(Ca0)2Si02 

C2S 

10-30 

Tri-Calcium  Aluminate 

CasA^Os 

(Ca0)3Al203 

C3A 

3-13 

Tetra-Calcium  Aluminoferrite 

Ca4Al2Fe20io 

(CaO)4.  Al203Fe203 

c4af 

5-15 

Calcium  Sulfate  Dehydrate 

(Gypsum) 

CaS04.2H20 

(Ca0)(S03)(H20)2 

c-s-h2 

3-7 

Source:  ref  [3], 


Cement  paste  is  made  by  combining  dry  Portland  cement  and  water.  The  process  of 
making  cement  paste  is  referred  to  as  hydration,  during  which,  there  is  loss  of  workability, 
solidification  and  hardening.  The  hydration  of  the  different  clinker  phases  results  in  different 
hydration  product  phases,  which  makes  up  the  cement  paste.  The  hydrated  cement  phases 
include;  calcium  silicate  hydrate  (C-S-H),  calcium  aluminum  hydrate  (C-A-H)  and  calcium 
hydroxide  (C-H)  [3,  4],  see  Table  2  for  the  common  components  of  hydrated  cement  paste. 
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Table  2 


Hydrated  Cement  Constituent 


Name  or  Mineral  Phase 

Chemical  Formula 

Cement  Chemical 

Notation  (CCN) 

Calcium  Silicate  Hydrate 

2(CaO).SiO2.0.9-1.25(H2O)  and/or 

CaO.Si02. 1 . 1(H20)  and/or 

0.8- 1 ,25(Ca0).Si02. 1 ,0-2.5(H2O) 

C-S-H 

Calcium  Hydroxide 

Ca(OH)2  OR  CaOH20 

CH 

Calcium  Aluminium  Hydrate 

More  complex  than  C-S-H 

C-A-H 

Aluminate  Ferrite  Trisulfate 

C3AS3H30-32 

Aft 

Aluminate  Ferrite 

C2ASHi2 

AFm 

Hydrogamet 

3Ca0.Al203.6(H20) 

c3ah6 

The  C-S-H  phase  which  constitutes  a  significant  percentage  of  the  hydrated  cement  paste 
is  of  interest  because  it  is  primarily  responsible  for  the  strength  and  load  bearing  attributes  of 
cement  [3,  5],  However,  because  of  its  complex  material  chemistry  [6],  C-S-H  molecular 
structural  representation  is  still  inconclusive,  and  hence  has  been  described  by  different 
molecular  model  crystal  structures  [7],  Some  of  these  molecular  structures  are  discussed  in 
chapter  3. 

The  interest  and  motivation  of  the  present  work  is  to  understand  the  effect  of  material  ion 
exchange  and  material  chemistry  level  changes  on  mechanical  stiffness  properties  and  shear 
deformation  behavior  obtained  through  computational  material  chemistry  level  modeling  of 
hydrated  cement  constituent  C-S-H  Jennite. 
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1.2  Cement:  Multi-Scale  Nature 

A  phase  in  continuum  micromechanics,  is  a  material  domain  that  can  be  identified  at  a 
given  length  scale  with  homogenous  deformation  state  and  constant  material  properties. 
Continuum  micromechanics  have  made  it  possible  to  represent  heterogeneous  regions  with 
homogeneous  phases,  with  equivalent  mechanical  properties  [8].  Individual  homogenous  phases 
at  microscopic  state  have,  on-average,  a  constant  strain  state;  and  classical  homogenization 
techniques  helps  to  predict  the  overall  behavior  of  heterogeneous  materials  from  that  of  their 
constituents  [9];  hence  the  possibility  of  homogenizing  the  lower  scale  into  higher  scale  and 
delivering  specific  scale  phase  properties,  volume  fraction  and  specific  morphologies  [3]. 

Cement  paste  exhibit  multi-scale  features;  it  has  a  complex  nature  with  random  microstructure  at 
different  length  scales  and  the  concept  of  the  continuum  micromechanics  have  been  used  in 
cement  structure  studies  [8-10].  Several  application  properties,  have  utilized  the  continuum 
micromechanics  length  scale  framework;  including  heterogeneous  thermoelasticity,  rate- 
independent  elastoplasticity,  nonlinear  elasticity,  viscoplasticity,  viscoelastic  coupling  and  rate- 
dependent  effects  [4],  these  are  primarily  at  the  microstructure  level.  However,  in  cement  paste, 
four  levels  of  scales  have  been  identified  [8]  i.e.  molecular  /nano-scale  level  determined  by  the 
material  chemistry  features,  micro-scale  with  microstructural  morphology  features,  meso-scale 
and  macro-scale.  Distinct  length  scales  and  features  are  based  on  the  requirement  that  each  scale 
should  be  separated  from  the  next  scale  by  at  least  one  order  of  length  magnitude,  which  is  also 
stated  in  the  literature  as  a  prerequisite  for  the  application  of  continuum  micromechanics  [3,4], 
As  a  result  of  the  multi-scale  nature  of  cement  paste,  the  properties  and  material  composition  of 
the  macro  and  engineering  scale  cement  can  be  controlled  or  affected  by  the  molecular  level.  A 
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potential  ability  to  alter  and  manipulate  the  chemistry  level  structure  of  hydrated  component  of 
traditional  C-S-H,  a  structure  predominantly  influenced  by  calcium  -  a  major  constituent  of  the 
cement  paste  at  the  fundamental  material  chemistry  level  -  will  imply  that  the  elastic  and  shear 
deformation  properties  of  the  paste  would  be  impacted.  Figures  1  and  2  shows  the  different 
levels  of  scale  associated  with  cementitious  materials.  Figure  1  shows  a  schematic  of  the  muti- 
scale  material  structure  while  Figure  2  shows  the  associated  forms  of  cement  in  engineering 
applications. 


Figure  1.  Schematics  of  multi-scale  material  structure  of  cement-based  materials, 
source  ref  [8]. 


Cement  Cement  Paste  Mortar  Concrete 


Figure  2.  Multiscale  nature  of  cement-based  materials. 
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1.3  Objectives 

Traditional  cement  paste  is  primarily  calcium-based,  but  engineering  community,  the  cement- 
producing  industry,  and  government  establishments  such  as  the  U.S.  Department  of  Defense  are 
looking  for  alternate  formulations  of  cement;  which  may  provide  tailored  properties  for  specific 
applications,  could  utilize  local,  and  cheaper  materials.  Another  factor  influencing  the  alternate 
forms  cement  cited  in  the  literature  is  the  need  to  reduce  the  energy  required  to  produce  the 
clinker  phase,  and  the  carbon  foot  print  of  the  cement  industry  [2],  Hence  there  have  been  efforts 
to  look  for  alternate  cement  formulations  that  are  not  calcium-based,  might  be  cheaper,  and  have 
the  potential  for  reducing  carbon  footprint  in  the  manufacturing  stage  of  calcium  based  cement. 
Several  forms  of  cement  formulations  has  been  sought,  for  instance,  other  forms  of  belitic 
cements  are  sought  by  trial  and  error  to  improve  the  reactivity  of  the  clinker  phase.  Some 
additional  approaches  that  have  been  used  trying  to  improve  the  properties  of  cement-based 
materials  include  thermal  processing  [11,  12]  and  addition  of  new  chemical  compounds  and  ions 
[2,  11,  12],  Most  common  types  of  ionic  substitution  which  has  been  used  in  cement  phases 
include;  Mg2+  for  Ca2+,  2A13+  and  2Fe3+  for  3Ca2+,  and  2Ca2+  for  Si4+  [2,  13].  The  present  work 
investigates  the  influence  of  replacement  of  calcium  ion  with  magnesium  ion  in  the  material 
chemistry  level  of  C-S-H  with  a  focus  on  mechanical  stiffness  and  shear  deformation 
characteristics  via  computational  material  modeling. 

Several  studies  have  focused  on  different  aspects  of  cement  material  improvement.  Hegoi 
and  co-workers  [2]  studied  the  effects  of  the  presence  of  chemical  substitution  on  the 

2_|_ 

physicochemical  properties  of  cement  clinker  phases  -  alite  and  belite.  They  incorporated  Mg  , 

O  I  O  I 

A1  and  Fe  into  the  structure  using  classical  forcefield  methods  and  reported  that  the 
crystallographic  site  within  the  unit  cell  is  equally  probable  for  Mg,  A1  and  Fe  substitution.  It  has 
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also  been  reported  that  Mg  incorporation  does  not  change  the  electronic  structure  appreciably 
compared  to  A1  and  Fe  incorporation  [2],  and  that  there  are  no  preferential  substitution  for  any 
calcium  sites,  when  Mg2+,  Al3+  and  Fe3+  were  incorporated.  The  main  focus  of  their  research  was 
on  the  reactivity  of  the  clinker  phase  for  reduced  energy  production.  Research  efforts  to  study  the 
effect  of  Mg  ion  exchange  on  the  mechanical  stiffness  and  deformation  have  not  been  reported. 

As  discussed  in  the  literature,  Mg  provides  a  feasible  incorporation  option  in  cement 
paste.  In  their  study  of  the  hydration  process  of  cement,  Stephan  and  co-workers  [13],  provided 
analytical  evidences  that  in  low  concentrations,  Mg  does  not  change  the  hydration  process  of 
C3S  -  a  component  of  cement  ( Table  1 ),  although  the  reactivity  may  increase  with  weight 
percent  and  with  cement  age.  Also,  because  of  the  relative  availability  and  the  ionic  similarities 
to  Calcium,  Mg  can  be  considered  as  a  potential  candidate  for  ionic  replacement.  Additionally, 
with  the  replacement  of  Ca  by  Mg  a  lower  temperature  for  processing  the  clinker  phase  can  be 
achieved. 

The  objective  of  this  work  is  to  understand  the  influence  of  the  effects  of  the  exchange  of 
Calcium  ion  with  Magnesium  ions  in  hydrated  cement  material  chemistry  structure  based  on  the 
shear  deformation  behavior  and  mechanical  stiffness  properties  of  the  modified  material  using 
computational  material  modeling  and  simulations.  Different  computational  simulation  methods 
have  been  employed  in  the  past  to  study  heterogeneous  materials.  Hegoi  et  al  [2]  used  a 
combination  of  forcefield  and  DFT  atomistic  simulation  in  their  mentioned  prior  mentioned 
study,  others  have  used  ab  initio  methods  to  study  the  structure  of  C-S-H  [14,  15].  There  are  also 
studies  of  the  mechanical  behavior  of  heterogeneous  materials  using  multi-scale  modeling,  and 
studies  to  understand  cement  hydration  with  finite  element  modeling  as  well  as  temperature 
dependency  of  the  microstructure  of  cement  hydrates  [8,  16-18]  reported  in  the  literature.  To  the 
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best  knowledge  of  the  author,  no  other  studies  have  been  focused  on  the  effects  of  chemistry 
level  ion  exchange  on  the  mechanical  stiffness  and  deformation  behavior  of  cement  hydrates. 

The  present  study  focuses  on  the  effect  of  Calcium  ion  substitution  with  Magnesium  on 
chemistry  level  structure  of  C-S-H  and  to  understand  the  variations  on  the  predicted  mechanical 
stiffness  and  shear  deformation. 

1.4  Computational  Modeling 

In  order  to  move  nanotechnology  research  from  the  nascent  exploratory  research  which 
was  initiated  mostly  to  understand  underlining  principles  to  designing  new  materials  that  will 
meet  existing  and  new  needs,  and  thus  become  commercially  viable;  a  rational  nanomaterial 
design  and  rigorous  engineering  design  would  be  required.  The  knowledge  and  understanding  of 
the  properties  based  on  atoms  and  molecules  that  form  the  material  chemistry  structure  at  a  nano¬ 
scale  is  very  critical  to  achieving  this  rational  design.  At  the  bottom  line  of  this  understanding  is 
the  use  of  computational  algorithms  and  relevant  analysis  code  developments  that  would  help  in 
predicting  the  behavior  of  new  material  performances  even  before  they  are  formed  and  subjected 
to  any  laboratory  efforts.  Material  chemistry  level  modeling  following  the  principles  and 
techniques  commonly  grouped  under  Computational  Material  Science  is  one  of  the  key 
facilitator  of  material  science  that  is  experiencing  fast  growth  pace  and  forms  the  basis  of  the 
present  work  [19]. 

The  background  of  computational  modeling  methods  and  framework  followed  in  the  present 
research  is  presented  and  discussed  in  Chapter  2. 
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CHAPTER  2 

Computational  Modeling  Methods  -  Background 

This  chapter  presents  the  background  details  of  the  computational  modeling  methods 
employed  in  the  present  work.  In  particular,  material  chemistry  level  modeling  is  based  on 
molecular  configurations  formed  with  atoms  and  bonds  and  their  interactions.  A  brief  discussion 
of  the  background  of  Molecular  Dynamics  (MD)  modeling  techniques  is  presented  next. 

2.1  Molecular  Dynamics  Techniques 

Quantum  mechanics  offers  an  accurate  solution  to  the  position-space  wave  function  of 
electron  level  systems.  However,  it  is  computationally  complex  and  difficult  to  solve  the  wave 
function  for  larger  systems  [20].  Hence,  Molecular  Dynamics  (MD)  methods  has  been  developed 
and  used  to  solve  atomic/  molecular  level  systems.  Unlike  quantum  mechanics,  molecular 
dynamics  can  be  applied  to  larger  and  complex  systems;  it  is  also  computationally  less  expensive 
and  faster.  Molecular  Dynamics  helps  to  fill  the  gap  created  when  statistical  mechanics  -which 
uses  partition  function  of  systems  in  equilibrium,  meets  a  roadblock  due  to  essentially  lack  of 
system  equilibrium.  This  is  because,  once  out  of  equilibrium,  theory  cannot  provide  quantitative 
solutions  unless  several  approximations  are  applied;  hence  simulations  like  Molecular  Dynamics 
are  employed  to  solve  this  equilibrium  challenges  [21].  A  particular  subset  of  Molecular 
dynamics  modeling  techniques  is  the  one  denoted  as  Classical  Molecular  Dynamics,  which  is  a 
numerical  technique  that  uses  classical  mechanics  to  model  atomic  scale  systems.  In  it,  the 
position  and  velocities  of  atoms  forming  a  molecular  system  are  estimated  based  on  the  solution 
of  the  classical  equations  of  motion,  Newton’s  laws  of  motion.  This  technique  cannot  be  used  to 
study  the  chemical  reactions  or  predict  chemical  reactivity  of  molecules. 
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Unlike  lattice-based  simulation  methods  like  cellular  automata,  MD  simulations  operates 
in  continuum  and  are  not  spatially  discrete  [22],  hence  problems  like  unfavorable  effects  due  to 
lattice  symmetry  and  lack  of  range  of  particle  velocities  are  not  encountered  in  MD,  although  it 
is  computationally  more  demanding  [21].  Other  simulations  techniques,  like  Monte-Carlo 
simulations,  eliminate  the  momentum  part  of  the  phase  space  and  considers  only  configuration 
space,  hence,  can  only  be  used  to  study  systems  in  equilibrium  and  process  consequences. 
However,  MD  enables  a  detailed  study  of  the  dynamics  (path  and  consequences)  of  a  molecular 
system  [23].  In  Classical  MD,  atoms  are  viewed  as  mass  points,  and  are  described  by  mass  ( m ), 
charge  (q),  position  (x,  y,  z)  and  velocities  ( Vx,  Vy,  Vz )  -  this  has  been  described  by  the  particle 
model  with  the  particles  carrying  the  properties  of  physical  objects  [20].  By  virtue  of  their 
existence  and  vicinity,  these  particles/atoms  interact  with  one  another  (bonds  can  be  viewed  as 
springs  and  the  bond  strength  are  related  to  the  spring  constant),  and  these  interactions  are 
referred  to  as  interatomic  potential.  Figure  3  shows  a  3-atom  system  with  their  various 
interactions.  In  molecular  dynamics,  the  interatomic  potentials  are  evaluated  via  force-field 
calculations  [21,  24], 


z 

A 
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Figure  3.  Atomic  positions  and  interaction  forces  on  a  3-atom  system. 
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2.2  Force  Field/Interatomic  Potential 

In  Classical  MD,  the  energy  associated  with  molecular  systems  is  based  upon  energy 
associated  with  the  atomic  level  interactions  of  different  atoms  in  a  molecular  system.  The 
potential  energy  of  atomic  system  is  expressed  as  the  sum  of  various  interactions  that  exist  in  a 
given  molecular  [25]  system  as; 

E Total  —  E Bonded  ~b  ^Non-bonded  ~b  E Cross  terms  Equation  1 

E Bonded  E stretch  ~b  E Angle  "b  ^Torsion  ~b  E  out— of  -plane  Equation  2 

Enon-bonded  Eyow  ~b  E Electrostatic  Equation  3 

where  ETotai  is  the  total  energy  of  an  atomic  system  accounted  for  by  a  combination  of  both 
bonded  and  non-bonded  components  of  the  interactions  of  the  atoms.  The  bonded  component  of 
the  interaction  potential  comprise  of;  the  bond  stretching  component  ( Estretch ),  which  views  the 
bonds  as  springs  under  harmonic  oscillation  and  thus  the  relation  in  figure  4  below.  The  bond 
angle  bending,  and  dihedral  angle  torsion  components  of  the  angular  distortion  ( EAngie  and 
ET or sion )  [25],  relate  the  interactions  between  three  and  four  atoms,  with  respect  to  the  angles 
between  three  atoms  and  the  angles  between  different  planes  of  four  atoms,  and  the  inversion 
term  ( E0ut_0f_piane ).  The  non-bonded  component  is  comprised  of  Vander  Waals  interactions 
(. EVDw )  which  relates  the  long  range  attractive  and  short  range  repulsive  effects  on  neutral  atoms 
using  the  Lennard-Jones  potential  [21,  25],  figure  4,  and  electrostatic  interactions  (EElectrostatic) 


[26], 
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Figure  4.  Bonded  and  unbonded  interactions. 

The  interatomic  force  can  then  be  obtained  as  the  spatial  gradient  of  the  total  potential  energy  of 
the  system  [20,  21,  27]. 

F(  =  —  ^  Equation  4 

The  trajectories  of  individual  particles,  rt(t)  of  the  system  can  be  estimated  solving  the 
Newton’s  equation  of  motion  [20,  21,  27]. 

d,ip  ■ 

Fj  =  mjCi;  =  rni~j^  =  mi  Equation  5 

During  the  time-evolution  of  the  dynamics  of  a  molecular  system  the  positions  and 
velocities  are  updated  after  each  time  step.  Since  all  the  thermodynamic  properties  of  the 
molecular  system  can  be  estimated  based  on  the  set  of  positions  and  velocities,  a  time-profile  of 
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them  can  be  created  during  the  time-evolution  of  the  molecular  system.  The  physical  properties 
under  thermodynamic  conditions  employed  in  the  simulations  can,  therefore,  be  average  of  large 
enough  periods  of  time  as  to  satisfy  the  conditions  required  by  the  Ergodic  hypothesis.  Figure  5 
shows  a  typical  molecular  dynamic  process  step. 
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Figure  5.  Molecular  dynamics  process. 

Figure  6  shows  a  time  average  line  for  a  thermodynamic  or  physical  quantity.  Section  2.5 
discusses  the  approaches  used  in  executing  this  transient  dynamic  process. 
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Figure  6.  Time  average.  Quantities  of  interest  are  obtained  as  average  over  a  dynamic  time. 
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The  solutions  to  the  equations  of  motion  are  based  on  constraints  such  that  certain 
thermodynamic  properties  are  controlled.  Hence,  the  system  under  investigation  can  be  subjected 
to  constraints  such  as  constant  temperature,  pressure  and  volume.  The  combination  of  these 
thermodynamic  constraints  is  called  ensembles  which  are  discussed  later.  In  Molecular 
dynamics,  the  potential  energy  for  bonded  and  non-bonded  interactions  of  different  material 
atoms  has  been  defined  by  different  forcefields  that  have  been  developed  and  applicable  to 
material  atoms  involved  in  the  molecular  system.  The  predicted  properties  of  a  molecular 
dynamics  analysis  for  material  modeling  are  influenced  by  the  underlying  forcefields  used  in 
analysis.  Two  such  force  fields  have  been  employed  in  the  present  study  on  material  modeling 
based  on  material  chemistry  structure  of  hydrated  C-S-H  Jennite. 

2.2.1  Universal  force  field  (UFF).  Universal  Force  Field  (UFF)  provides  the  possibility 
of  constructing  and  defining  most  structural  features  across  the  elements  in  the  periodic  table. 

The  limitation  of  prior  existing  forcefield  applications  to  proteins,  nucleic  acids  and  organic 
molecules  deepened  the  efforts  resulting  in  the  development  of  the  Universal  force  field  which 
could  be  applied  to  the  entire  periodic  table.  It  is  generated  using  parameters  which  include 
hybridization  dependent  atomic  radii,  Van  der  Waals  parameters,  hybridization  angles,  effective 
nuclear  charges,  as  well  as  inversion  and  torsional  barriers.  UFF  has  126  atom  types,  which  are 
described  by  a  five  character  mnemonics.  The  force  field  assigns  potential  energy  related  to  the 
geometry  of  the  molecules  as  a  superposition  of  various  two,  three  or  four  body  interactions.  For 
atoms  that  are  bonded  to  each  other  in  a  1,  2  interaction  or  those  bonded  to  a  common  atom  in  a 
1,  3  interactions,  UFF  follows  the  convention  of  excluding  Van  der  Waals  and  electrostatic 


interaction  [28]. 
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2.2.2  COMPASS  force  field  (COMPASS).  Condensed-phase  optimized  molecular 
potentials  for  atomistic  simulation  studies  or  COMPASS  force  field  was  developed  as  a  general 
all-atom  forcefield  for  atomistic  simulations  of  common  organic,  small  inorganic  molecules  and 
polymers.  This  is  a  proprietary  forcefield  within  the  MD  analysis  and  visualizer  code  Materials 
Studio.  Ab-initio  and  empirical  parameterization  techniques  were  used  [29],  and  validated  using 
condensed-phase  properties  for  molecules  in  isolation.  This  force  field  can  be  used  for  a  broad 
range  of  molecules  in  isolation  and  in  condensed  phases.  It  is  capable  of  predicting  numerous 
solid  state  properties,  unit  cell  structures,  lattice  energies,  elastic  constants  and  vibrational 
frequencies  [30]  and  under  a  wide  range  of  conditions  of  temperature  and  pressure. 

2.3  Ensembles  and  Controls 

Molecular  dynamics  is  based  on  the  integration  of  classical  equations  of  motion  on 
molecular  systems  to  simulate  their  kinetic  and  thermodynamic  properties.  It  samples  micro- 
canonical  ensembles  (NVE)  [31],  since  a  system  is  usually  characterized  by  a  time-independent, 
translationally  and  rotationally  invariant  Hamiltonian.  The  integration  of  the  classical  equations 
of  motion  leads  to  a  trajectory  that  maps  several  dynamic  micro-canonical  ensembles  (NVE), 
where  all  the  forces  which  appear  in  the  Newton’s  equations  of  motion,  are  related  to  the 
potential  energy  of  the  system,  and  the  total  energy  of  the  system  ( ETotai  =  EKinetic  + 

E potential )  is  conserved  (NVE  -  fixed  number  of  particles,  volume  and  energy).  But  there  are 
cases  when  it  is  desirable  to  make  a  system  comparable  with  experiments.  For  instance  when 
average  temperature  must  be  at  desirable  values,  this  leads  to  a  canonical  sample  configuration 
with  fixed  number  of  particles  and  volume  and  with  the  temperature  having  a  specified 
macroscopic  average  but  with  fluctuating  instantaneous  total  energy  or  Hamiltonian  (NVT 
ensembles).  Other  possibilities  include,  grand  canonical  ensembles  (/iVT)  where  there  is  a 
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constant  volume  and  temperature  as  in  canonical  ensembles  but  it  is  however  open  to  particle 
exchange  with  the  surrounding  bath.  In  grand  canonical  ensembles,  the  chemical  potential  of 
different  species  have  a  specified  average  [31].  In  isothermal-isobaric  (Gibbs)  ensembles  (NPT), 
the  pressure  has  a  specified  average  and  the  instantaneous  volume  of  the  system  can  fluctuate 
[31].  Thermostats  algorithms  are  the  various  computational  based  methods  in  molecular 
dynamics  which  are  introduced  to  modify  the  time  dependent  changes  in  atom  positions  and 
velocities  to  achieve  this  temperature  controls  [31,  32],  The  purpose  of  these  thermostat  controls 
is  to  ensure  the  average  temperature  of  a  molecular  system  is  maintained  at  and  close  to  a  desired 
temperature  [32],  A  brief  background  of  the  different  thermostats  employed  in  molecular 
dynamics  analysis  and  studied  in  the  present  work  is  presented  next. 

2.3.1  Andersen  thermostat.  It  works  by  coupling  the  molecular  system  to  a  heat  bath,  by 
an  occasional  stochastic  collision  that  acts  randomly  on  selected  particles,  imposing  a  desired 
temperature  on  the  system  [32],  The  instantaneous  stochastic  collision  affects  the  momentum 
(velocity)  of  one  particle  and  other  particles  are  unaffected  by  the  collision  [32,  33].  The 
Hamiltonian  equations  for  all  the  particles  are  integrated  until  the  time  for  next  instantaneous 
collision  [32-34],  Andersen  thermostat  can  be  used  for  non-equilibrium  processes  where  large 
amount  of  heat  is  liberated  for  systems  with  low  heat  capacity,  and  for  processes  such  as 
nucleation  and  phase  separations  which  require  large  energy  and  density  fluctuations  [34], 

2.3.2  Velocity  scaling.  It  is  probably  the  most  straight  forward  temperature  control 
method.  Maxwell-Boltzmann  distribution  is  used  to  draw  the  distribution  of  velocities  of 
material  components.  The  average  kinetic  energy  of  each  component  per  degree  of  freedom  can 
be  obtained  and  related  to  the  temperature  using  equipartition  theorem,  and  this  is  averaged  over 
all  the  particles  to  obtain  the  instantaneous  temperature  of  the  finite  sized  material  system.  If  the 
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temperature  used  to  generate  the  velocity  distribution  does  not  coincide  with  the  instantaneous 
temperature,  then  the  velocity  is  rescaled  until  the  temperature  coincides  [35].  Velocity  scaling 
does  not  reproduce  the  canonical  ensembles  [35,  36],  and  it  is  not  deterministic  or  time 
reversible  and  does  not  remove  localized  and  unwanted  correlation  motions  [36]. 

2.3.3  Nose  thermostat.  Instead  of  using  stochastic  collision  to  introduce  energy 
fluctuations  into  the  simulation  system,  a  Lagrangian  with  artificial  coordinates  and  velocities 
also  referred  to  as  extended  Langragian  method  are  used  in  Nose  thermostat.  A  Nose  thermostat 
controlled  system  has  more  independent  variables  than  equivalent  statistical  mechanics  approach 
giving  more  accurate  results  for  static  quantities  [32,  37,  38].  It  has  also  been  shown  that  the 
simulations  of  extended  systems  (which  includes  a  real  micro-canonical  system  and  a  heat  bath) 
produces  a  canonical  ensemble  in  the  real  systems  due  to  heat  exchange  between  fictitious 
degree  of  freedom  and  real  systems,  with  coupling  between  the  two  controlled  by  a  fictitious 
mass  [35].  The  method  offers  a  stable  and  efficient  approach  to  molecular  dynamics  analysis 
simulations.  However,  its  more  expensive  to  carry  out  each  optimization  steps  [32],  and  it  is  not 
Ergodic  (i.e.  the  system  does  not  have  has  a  finite  chance  to  go  everywhere  in  phase-space  in  a 
finite  time)  in  some  difficult  cases  such  as  in  harmonic  systems  [33]. 

2.3.4  Berendsen  thermostat.  Berendsen  et  al  [39]  introduced  weak  coupling  of  system 
to  heat  bath  thereby  limiting  the  temperature  fluctuations  in  canonical  ensembles  found  in 
velocity  scaling.  In  this  temperature  control  method,  the  coupling  adds  or  removes  energy  form 
the  system  to  maintain  the  temperature  and  the  velocities  are  scaled  at  each  step  such  that  the  rate 
of  temperature  change  is  proportional  to  the  difference  between  the  target  and  instantaneous 
temperature  given  by; 


Equation  6 
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When  the  target  temperature  To  is  less  than  desired  system  temperature  T,  the  temperature 
will  increase  and  if  T>T0,  heat  will  be  removed  from  the  system.  A  generally  employed  value  for 
x  is  0.4ps,  which  result  in  modest  temperature  fluctuations  in  most  MD  analysis  simulations  [35, 
39].  Berendsen  thermostat  leads  to  exponential  relaxation  of  the  temperature  to  a  target  one,  it 
does  not  strictly  fix  the  temperature  [35]  .  Also  when  the  coupling  is  weak  (r  =  O.Olps )  there 
may  be  perturbation  of  the  various  quantities  that  will  not  lead  to  correct  canonical  averages  [35]. 

2.3.5  Andersen  barostat.  The  methods  in  MD  by  which  the  time  dependent  changes  in 
atomic  configurations  can  be  modified  as  a  function  of  pressure  are  called  the  barostat.  Andersen 
barostat  mimics  the  action  of  a  piston  by  coupling  the  system  to  external  variable  [27],  by 
addition  of  an  extra  degree  of  freedom  to  the  volume  of  the  simulation  box  as  well  as  associated 
potential  and  kinetic  energy  terms.  It  allows  for  volume  change  within  the  simulation  with 
average  volume  determined  by  a  balance  between  the  internal  pressure  of  the  system  and  the 
desired  external  pressure  of  the  system. 

2.3.6  Parinelo  barostat.  Parinelo  barostat  controls  the  pressure  by  the  introduction  of 
additional  degrees  of  freedom  to  the  whole  system  of  atoms  in  time  and  space,  via  the 
transformation  of  the  spatial  coordinates  into  scaled  positional  coordinates,  x£  e  [0,1)3.  This  is 
possible  by  transforming  the  Lagrangian  describing  the  system  into  the  corresponding 
Hamiltonian  [20], 

2.3.7  Berendsen  barostat.  Similar  to  the  Berendsen  thermostat,  the  barostat  is  coupled 
with  a  pressure  bath  by  a  weak  coupling  with  coupling  constant  rp  and  the  volume  of  the  system 
is  scaled  by  a  scaling  factor  X.  The  rate  of  change  of  pressure  is  given  by  Equation  7  and  scaling 
factor  is  defined  by  Equation  8. 

=  p  ( Pbath  ~  P(t))  Equation  7 
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1 

1  =  1  —  k  —  =  —  (P  —  Pbath)  Equation  8 

Tp  Tp 

where;  Pbath  is  the  pressure  of  bath,  P(t)  is  the  actual  pressure  at  time  t  and  atomic  coordinates  are 
scaled  by  a  factor  of  A1^,  new  atom  position  is  defined  as; 

riscaie  =  ^3ri  Equation  9 

2.4  Energy  Minimization 

Consider  a  simulation  box  with  the  atoms  of  the  molecular  structure  located  on  certain 
defined  locations  in  the  x,y,z  space.  The  structures  that  define  the  molecular  material  model 
configuration  are  assigned  specific  locations  for  each  of  their  atoms  which  may  not  necessarily 
be  the  equilibrium  positions  with  a  minimal  energy.  In  order  to  evaluate  the  energy  of  the 
system,  we  need  a  force-field  or  interaction  potential  which  is  nothing  but  a  mathematical 
expression  of  the  relationships  that  are  dependent  upon  the  atom  type,  locations  of  each  atoms 
and  the  resulting  energy  based  on  bonded  and  non-bonded  interactions.  Hence,  the  potential 
energy  is  E  =  E(X()  in  general.  In  the  present  work,  we  use  the  UNIVERSAL  and  COMPASS 
force  fields  to  relate  the  above  interactions.  Before  any  dynamics  run  in  a  molecular  dynamics 
analysis,  energy  minimization  is  employed  to  get  a  low  energy  stable  configuration  of  the 
material  system.  There  are  numerous  mathematical  algorithms  that  can  be  employed  for  energy 
minimization  of  the  multi-dimensional  energy  function:  conjugated  gradient,  Newton-Rapson 
iteration  etc.  The  numerical  integration  method  used  for  solving  the  classical  equation  of  motion 
is  presented  next. 

2.5  Time  Integration  Method 

The  integration  of  the  equations  of  motion  is  done  numerically  and  there  are  different 
methods  of  doing  this  in  principle  [21,  40],  However,  there  are  selection  criteria  which  include 
limited  calculation  of  the  force  component,  a  corresponding  increase  in  time  step  if  force 
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components  are  calculated  but  without  crossing  the  upper  bound  limit  constraints  placed  by 
Lennard-Jones  repulsive  force  at  short  distances  [21].  Hence  Runge-Kutta  and  adaptive  methods 
of  time  integration  which  could  not  increase  their  time  steps  beyond  the  Lennard-Jones  limits  are 
quickly  discarded  in  most  MD  analysis.  Of  common  use  are  the  leap  frog  methods,  Verlet 
methods  and  the  predictor  corrector  methods. 

Leap  frog  and  Verlet  method  are  algebraically  equivalent  [21]  and  are  perhaps  the  most 
widely  used  time  integration  methods  in  molecular  dynamics,  easy  to  implement,  stable  and  time 
reversible.  Verlet  method  is  used  in  the  MD  simulations  performed  in  this  study.  It  is  a  direct 
solution  to  the  second  order  equation  5  [21,  27].  It  follows  from  the  Taylor’s  series  expansion  of 
coordinates  and  based  on  the  previous  particle  position  r(t  —  St),  present  position  r(t)  and 
acceleration  a(t).  The  next  position  can  be  written  without  necessarily  evaluating  velocity  the 
term  as: 

r(t  +  St)  =  2  r(t)  —  r(t  —  St)  +  8t2a(t)  Equation  10 

The  velocity  term  which  is  not  needed  for  trajectory  calculations  have  been  eliminated  by 
Taylor  series  expansion  about  r(t).  It  can  however  be  computed  when  needed  for  kinetic  energy 
calculations  as  [21,  27], 

V(t)  =  r(t+<5t^(t~<5t)  Equation  11 

2.6  Periodic  Boundary  Condition 

Finite  and  infinite  systems  behave  in  different  manner  [21].  The  aim  of  Periodic 
Boundary  conditions  (PBC)  in  the  molecular  dynamics  analysis  is  to  ensure  that  the  molecular 
system  does  not  have  an  abrupt  boundary  with  a  vacuum,  to  remove  any  significant  peripheral 
surface  effects  [30],  In  order  to  capture  the  typical  state  of  interior  atoms,  wall/boundary  effects 
are  eliminated  [21]  using  PBC.  Periodic  boundary  can  be  viewed  as  an  array  of  infinite  identical 
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copies  of  a  simulation  region  implemented  in  order  to  reduce  the  fraction  of  atoms  near  the 
walls/boundary  [27,  41].  Particles  that  leave  a  simulation  boundary  reenters  through  the  opposite 
face,  hence  the  boundary  of  the  central  cell  are  not  rigid  and  there  are  no  particle  reflection  in  the 
central  cell  [20],  but  the  periodic  images  of  the  central  cell  behaves  in  the  same  way.  It  is 
important  to  know  if  the  properties  of  both  macroscopic  and  small  infinite  periodic  system  are 
the  same.  The  common  practice  in  computational  simulation  is  to  ensure  the  choice  of  periodic 
size  has  little  effect  on  the  equilibrium  thermodynamic  properties  of  the  structure,  and  depending 
on  computational  resources  available,  to  increase  the  box  size  to  maintain  constant  density. 

Cubic  box  is  usually  used  in  most  periodic  repetition  [27]. 

2.7  Dynamic  Simulations  and  Post- Analysis  of  Physical  Quantities 

The  positions  and  velocities  of  molecular  system,  as  a  result  of  their  molecular 
interactions,  carry  information  that  can  be  used  to  evaluate  their  instantaneous  thermodynamic 
properties  such  as  their  temperature,  pressure,  kinetic  energy  and  potential  energy  (Equations  1-5 
and  figure  4).  From  section  2.5  above,  we  saw  that  we  can  evaluate  new  positions  of  particles 
based  on  the  numerical  time  integration  methods.  Hence,  at  each  time  steps  the  required 
thermodynamic  properties  can  be  computed  and  profiled. 

The  positions  and  velocities  of  the  molecular  system  in  their  initial  state  are  chosen  such 
that  no  particles  occupy  exactly  the  same  position  [21,  27,  42],  Lattice  locations  are  often  used  to 
place  particles  and  velocity  components  are  attributed  to  the  particles  [42]  from  Gaussian 
distribution  (Equation  13),  with  magnitudes  conforming  to  the  required  temperature  of  the 
system.  This  is  then  corrected  such  that  overall  momentum  is  zero  [27]. 

P  =  Yli=imiVi  =  0  Equation  12 


Equation  13 
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where  p(yix)  is  the  probability  density  function  of  the  velocity  component  vix  [27], 

Most  times  the  simulation  is  started  from  disordered  configuration  at  different  density  and 
temperature.  Equilibration  is  done  to  bring  the  molecular  system  to  a  new  equilibrium  state.  This 
is  done  up  to  a  satisfactory  time  where  there  is  small  oscillation  of  tracked  properties  such  as 
potential  energy  about  a  steady  mean  value  [27]. 

The  interatomic  potential  is  dependent  on  the  positions  and  velocities  of  the  atoms.  As  shown  in 
the  equations  in  figure  4. 

Also,  the  kinetic  energy  of  the  molecular  system  can  be  derived  from  equation  14  [27]; 

K.E  =  Equation  14 

Base  on  the  velocities  of  the  atoms,  the  molecular  system  temperature  can  be  obtained 
from  the  average  of  a  kinetic  temperature  function  derived  from  Virial  theorem  and  equipartition 
principle  for  N  atoms  with  Nc  internal  molecular  constraints  as  [27,  43]; 


T  = 


2-k 


(3  N-Nc)kB  (3N-Nc)kB 


hzuz2i=i\Pt\2/mt 


Equation  15 


where;  A  is  the  kinetic  energy  in  phase  function,  kB  is  the  Boltzmann  constant  and  Pj  is  the 
momentum  of  the  particle.  The  pressure  can  also  be  obtained  from  the  average  of  the  pressure 
function  [27,  43,  44], 

p  =  pkBT  +  W /V  Equation  16 

where;  W  is  defined  as  the  internal  virial  restricted  for  intermolecular  interactions,  and  V  is  the 
potential  energy  in  phase  function. 

Virial  expression  can  also  be  used  to  obtain  the  internal  stress  and  strain  tensors  of 
molecular  systems  in  atomistic  simulations  [30,  45],  these  are  also  based  on  the  position  and 
velocities  of  the  atoms  ( Equations  17,  18). 
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o-  =  "  ^  [(Sf  m,  (Vf V?-))  +  Ci<y  n,;  fi.j)  ]  Equation  1 7 

where  V0  is  the  un-deformed  system  volume  and  fi  is  the  force  acting  on  particle  i.  In  order  to 
obtain  the  stress  tensor  of  the  molecular  structure  in  its  static  state,  applied  force  stress  term, 

(2i<  /  ri,j  fi,j)  in  Equation  17  is  omitted  [30,  45].  The  virial  expression  for  the  strain  tensor  is 
given  by; 

£  =  ~  [(/io)_1(hrh)/to 1  —  1]  Equation  18 

where  h0,  and  h  are  the  matrices  formed  by  the  column  vectors  defining  the  material  system  at  a 
reference  state  and  at  corresponding  new  state.  When  Hooke’s  law  is  applied  on  the  stress  and 
strain  tensors,  elastic  compliance  and  stiffness  matrices  can  be  obtained  for  the  molecular 
material  configuration  [30]. 

Cijki  =  Equation  19 

derivations  of  how  the  elastic  constants  of  continuum  particle  can  be  obtained  from  the  material 
strain  tensor  for  either  isothermal  or  adiabatic  process  can  be  obtained  in  dedicated  textbooks 
[45].  The  next  section  gives  the  summary  of  the  elastic  constraint  matrix  used  for  elastic  stiffness 
property  determinations. 

2.8  Predictive  Mechanical  Properties 

Computational  simulation  enables  prediction  of  the  mechanical  stiffness  properties  of 
materials  (bulk  modulus,  shear  modulus  and  young’s  modulus)  using  their  elastic  constants 
which  can  be  determined  using  [45,  46]; 

Cjj  =  -  -—7—  Equation  20 

lJ  V  dSidSj  ^ 

Where:  £b  sj  are  lattice  strain  components,  U  is  potential  energy,  V  is  simulation  cell  volume. 

The  6><6  elastic  stiffness  matrix,  C  can  be  obtained  from  the  above  relation. 
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And  elastic  compliance  matrix,  S  =  C’1  can  also  be  obtained. 


Equation  21 


Equation  22 


W>61  ‘-*66' 

Hence  the  mechanical  stiffness  properties  for  the  homogeneous  stiffness  modulus  can  be 
predicted  following  Voight,  Reuss  or  Hill  elastic  stiffness  relations  from  the  elastic  and 
compliance  matrix  coefficients  as  [47] : 


_  1 

R  (^ll+^22+^33)  +  2(^12+^23+^3l) 

Equation  23 

Tyr  (C11  +  C22+C33)  +  2(C12+C23  +  C31) 

J\v  — 

V  9 

Equation  24 

r,  (KV+Kr) 

Kit  — 

M  2 

Equation  25 

Q  _  ^ 

R  4(S11+S22+S33)—4(S12+S23+S31)  +  3(S4.4.+S3$+S66) 

Equation  26 

r  (Cii+C22  +  C33)-(^12+^23+^3l)  +  3(Q4 +  ^55+^66) 

v  15 

Equation  27 

r  (gV  +  Gr) 

Uu  — 

M  2 

Equation  28 

3K-2G 

V  =  - 

6K+2G 

Equation  29 

E  =  2G(l  +  u) 

Equation  30 

where;  Kr=  Reuss  bulk  modulus,  Ky  =  Voight  bulk  modulus,  Kh 

=  Hill  bulk  modulus, 

Gr=  Reuss  shear  modulus,  Gv  =  Voight  shear  modulus,  GH  =  Hill  shear  modulus,  E=  Young’s 
Modulus,  v  =  Poisson  ratio  [48]  and  Cy  and  Sy  are  the  components  of  the  stiffness  and 
compliance  matrix  respectively  [47,  49] . 
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2.9  Estimation  of  Material  Property  Based  on  Deformation  and  Stress-Strain  Behavior 

Most  of  the  present  work  in  molecular  dynamics  modeling  of  mechanical  properties  is 
determination  of  stiffness  modulus  as  discussed  earlier.  Material  characteristics  are  defined  by 
their  stress-strain  behavior  under  deformation  loading  and  the  elastic  properties  of  a  material 
during  deformation  can  be  obtained  from  the  stress-strain  curve.  Subjecting  a  material  to  shear, 
compressional  or  tensile  force  are  ways  of  implementing  such  deformation.  For  material  under  a 
loading  condition,  an  increase  in  the  load  results  in  an  increase  of  the  strain,  first  linearly 
following  Hooke’s  law  in  the  elastic  deformation  region,  after  which  a  slippage,  non-linear, 
plastic  deformation  occurs  [48]  that  causes  response  that  is  observed  in  a  typical  stress-strain 
curve,  as  shown  in  figure  7. 


Figure  7.  Typical  stress-strain  response 
Source:  Ref  [48] 

Computational  material  modeling  of  this  stress-strain  behavior  based  on  material 
chemistry  is  still  limited  and  has  not  been  fully  studied  and  understood.  Recently,  shear 
deformation  versus  strain  behavior  was  recently  investigated  in  our  research  group  [50].  Other 
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recent  literature  had  studied  the  unidirectional  tensile,  quasi-static  shear  deformation  for 
nanoscale  C-S-H  Tobermorite.  Molecular  Dynamics  modeling  of  the  stress-strain  behavior  based 
on  material  chemistry  molecular  structures  provides  a  fundamental  methodology  to  understand 
the  predicted  behavior  of  a  material  solely  based  on  its  material  chemistry,  and  provides  an 
avenue  to  investigate  the  effect  of  the  changes  in  the  material  chemistry  on  the  mechanical 
behavior. 

Following  the  molecular  dynamics  modeling  methodology  established  in  our  research 
group,  the  effect  of  material  chemistry  of  Magnesium-modified  Jennite  on  shear  deformation  of 
nanoscale  traditional  Jennite  is  investigated  in  the  present  study.  In  the  current  work,  shear 
deformation  is  implemented  by  deforming  the  molecular  structure  by  increasing  a  particular 
angle  of  the  crystalline  structure  while  adjusting  the  length  of  the  system  to  guarantee  that  the 
volume  remain  constant  during  the  deformation  process.  Figure  8  shows  a  triclinic  structure 
under  shear,  along  the  c/i-plane  by  a  gradual  change  in  the  angle- a,  thus  subjecting  the  structure 
to  deformation  with  resulting  shear  strain  given  by  ecb  =  8acb.  For  every  increment  of  the 
angle- a,  the  material  experiences  a  shear  strain ,  8a  given  by  Equation  31. 


m  =  <Xi-i  +  5a 


Equation  31 
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Figure  8.  Shear  deformation 

The  resulting  stress  strain  curve  from  the  computation  model  can  then  be  used  to  estimate  and 
predict  the  stress-strain  behavior  and  elastic  properties  of  the  material  under  deformation  loading 
for  all  the  three  crystallographic  planes  as  discussed  later  in  chapter  4. 
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CHAPTER  3 

Computational  Material  Modeling  Analysis  of  C-S-H  Jennite 

3.1  Software  Analysis  Codes  for  Molecular  Dynamics  Modeling 

Different  molecular  dynamics  modeling  analysis  codes  have  been  developed  and  are 
available  [21].  Some  of  which  are  open  source  codes  or  commercial  codes.  Large-scale 
Atomic/Molecular  Massively  Parallel  Simulator  (LAMMPS)  and  GROningen  MAchine  for 
Chemical  Simulations  (GROMACS)  are  examples  of  open  source  MD  modeling  analysis  codes. 
Accelrys  Material  Studio  is  an  example  of  commercial  analysis  code  [21].  In  the  present  work, 
MD  modeling  computational  analysis  was  conducted  using  Materials  Studio. 

3.2  C-S-H  Gel 

As  mentioned  in  chapter  1,  C-S-H  gel  is  a  major  component  in  the  cement  paste  that  is 
responsible  for  the  load  bearing  attributes  and  strength  of  cement.  C-S-H  has  a  complex  material 
chemistry  [6].  The  molecular  structural  representations  of  the  C-S-H  gel  are  still  inconclusive, 
hence  more  than  one  structural  representations  are  available  [7].  Some  of  these  molecular 
representations  include: 

•  Wollastonite  group  which  comprises  of  Foshagite  (Ca^SisCLXOH^)  [51], 
Hillebrandite  (Ca2(Si03)(0H)2)  [52],  Xonotlite  (Ca6Si60i7(0H)2)  [53],  Okenite 
([Ca8(Si60i6)(Si60i5)2(H20)6]4'[  Ca2(H20)9)  3H20])  [54],  and  others  that  are 
compiled  by  Richardson  [7] 

•  Tobermorite  group  which  comprises  of  Clinotobermorite0  (CasSieOn  5H20), 
Clinotobermorited  (CasSieOn  5H20),  Clinotobermorite  9  A0,c  (CasS^O^  (OH)2), 
Clinotobermorite  9  A0,d  (CasS^O^  (OH)2)  [55,  56],  anomalous  and  normal 
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Tobermorite  11  A0  (Ca4Si60i5  (OH)2  5H20  and  Ca4.5Si60i6(0H)25H20 
respectively)  [57,  58]  and  Tobermorite  14A°  [59] 

•  Jennite  group  comprising  Jennite  and  Metajennite  (Ca9Si60ig(0H)6  8H20).  Most 
of  the  crystal  structures  are  similar  to  one  another  than  might  seem  apparent  at 
first  look,  and  a  complete  composition  and  crystal  information  for  most  of  the  C- 
S-H  and  other  related  phases  has  been  tabulated  by  Richardson  [7], 

From  the  various  molecular  chemistry  structural  representations,  Tobermorite  14 A0  and 
Jennite  structures  are  the  traditionally  accepted,  widely  used  and  adapted  molecular  structural 
representation  for  most  C-S-H  cement  paste  studies  [4,  6,  7,  60,  61].  For  this  study,  Jennite 
structure  of  the  material  chemistry  configuration  of  hydrated  cement  C-S-H  was  used. 

3.2.1  Jennite.  Jennite  is  a  representative  mineral  form  of  calcium  silicate  hydrate  (C-S- 
H).  It  has  a  chemical  formula  of  Ca9Si6018(0H)e.  8(H20)  [7,  60,  61].  It  was  named  after 
Clarence  Marvin  Jenni  (1896-1973)  for  its  discovery  [62],  Jennite  is  composed  of  single  chains 
with  repeating  unit  of  three  Si04  (Wollastonite-type  dreier  single  chains),  ribbons  of  edge 
sharing  CaC>6  octahedral  and  additional  CaC>6  octahedral  on  inversion  centers  [60,  62],  Jennite 
has  a  triclinic  structure  {figure  9)  with  dimensions  of  a  single  unit  cell  shown  in  table  3  [60].  The 
MD  modeling  analysis  employed  in  the  present  study  is  based  on  the  material  chemistry  structure 
configuration  of  C-S-H  Jennite. 

Table  3 

Cell  Parameters  of  a  Single  Unit  Cell  Jennite 


a 

b 

c 

a 

P 

y 

10.576A° 

7.265AU 

10.931A° 

101.300° 

96.980° 

109.650° 
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Key: 

Green  =  Ca 


White  =  H 
Gold  =  Si 
Red  =  O 


Figure  9.  Jennite  structure 

The  MD  modeling  analysis  methodology  employed  in  the  present  study  with  Accelrys  Materials 
Studio  is  briefly  presented  next. 

3.3  Molecular  Dynamics  Modeling  Process 

To  initiate  a  MD  simulation,  the  initial  positions  of  the  atoms  forming  the  molecular 
structure  is  subjected  to  energy  minimization  and  geometric  optimization  using  Material  Studio. 
The  minimal  energy  configuration  is  then  subjected  to  thermodynamic  pressure  and  temperature 
constrains  followed  by  a  dynamic  analysis  (time  evolution  of  the  system).  The  flow  diagram  of 
the  general  MD  simulation  scheme  is  presented  in  figure  10. 
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Figure  10.  MD  material  modeling  process 

The  present  study  focuses  on  the  understanding  the  effect  of  material  chemistry  changes 
in  C-S-H  Jennite  via  MD  analysis  and  development  of  relevant  MD  based  material  modeling 
methodology.  The  exchange  of  Calcium  ions  with  Magnesium  as  discussed  earlier  was  achieved 
in  MD  modeling  with  Accelrys  Materials  Studio  as  follows. 

3.3.1  Molecular  Structural  Magnesium  Exchange  Approach.  The  traditional  Jennite 
unit  cell  molecular  structure  was  obtained  from  the  American  Mineralogist  Crystal  Structure 
Database  [63],  This  initial  molecular  structure  configuration  was  loaded  into  Material  Studio  and 
was  labeled  as  JO.  For  supercell  structures  to  be  made  in  Material  Studio,  a  PI  space  group  must 
be  used  to  impose  the  symmetry  and  was  followed.  The  C-S-H  Jennite  unit  cell  structure 
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contains  9  calcium  ions.  The  Magnesium-modified  C-S-H  Jennite  molecular  structure  was 
obtained  by  successively  replacing  the  calcium  atoms  by  magnesium  atoms.  The  Magnesium- 
modified  molecular  structures  were  then  employed  in  the  present  work  to  understand  the 
influence  of  the  ion  exchange. 

A  systematic  study  of  replacing  the  different  number  of  Magnesium  ions  in  place  of 
Calcium  ions  was  conducted  in  the  present  work  starting  with  the  replacement  of  2  Magnesium 
ions.  Manzano  et  al,  had  established  that  there  are  no  significant  preferential  tendencies  about 
where  the  Ca  ions  are  replaced  by  Mg  ions,  i.e.  a  random  substitution  pattern  in  the  MD  model 
accounts  for  the  molecular  structures  expected  in  physical  experimentation  [2],  In  the  first 
modified  C-S-H  Jennite  configuration  with  Magnesium  atoms,  two  calcium  ions  were  randomly 
replaced  with  magnesium  ions,  creating  a  modified  Jennite  structure,  Jl.  The  two  Magnesium 
ions  configuration  resulted  in  a  14.8%  Mg/(Mg+Ca)  ratio.  Additional  Magnesium  modified 
structures  were  created  by  the  exchange  of  4  -  9  calcium  ions  resulting  in  different  magnesium 
concentrations  by  weight.  The  various  Magnesium  modified  molecular  structural  configurations 
of  C-S-H  Jennite  considered  for  the  present  study  are  shown  in  table  4. 

Table  4 

Traditional  and  Modified  Jennite  Structures  (%  Weight) 


Jennite 

n-Ca 

n-Mg 

n-H 

n-O 

n-Si 

Atomic  Mass 

“/.Weight  =  ^ioo 

JO 

9 

0 

20 

33 

6 

1077.35 

0.0 

Jl 

7 

2 

20 

33 

6 

1045.80 

14.8 

J2 

5 

4 

20 

33 

6 

1014.25 

32.7 

35 


Table  4 
Cont. 


J3 

3 

6 

20 

33 

6 

982.71 

54.8 

J4 

1 

8 

20 

33 

6 

951.16 

82.9 

J5 

0 

9 

20 

33 

6 

935.39 

100.0 

Molecular  Mass:  Ca  =  40.08,  Mg  =  24.31,  H  =  1.01,  0  =  16.00,  Si  =  28.09 

Figure  11  shows  the  traditional  and  modified  structures  of  the  single  unit  cell  of  C-S-H  Jennite 
with  different  color  representations  to  distinguish  the  atoms. 


Traditional  .Tennite  Modified  Jennite  (14.8%M 


Modified  Jennite  (54.8%Mg)  Modified  Jennite  (100%Mg) 


Green=Ca,  Blue=Mg,  White=H,  Gold=Si,  Red=0 


Figure  11.  Traditional  and  modified  Jennite  structures. 
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A  systematic  MD  analysis  study  of  these  different  molecular  structures  was  conducted.  As 
discussed  earlier,  several  computational  MD  simulation  parameters  were  determined  for  this 
material  system  prior  to  full  MD  analysis.  These  different  MD  parameters  that  were  investigated 
and  established  is  discussed  next. 

3.4  Parameter  Determination  Analysis 

There  are  several  parameters  that  need  to  be  established  for  specific  material  system 
when  using  MD  simulations.  Further,  in  Material  Studio,  the  configurations  for  the  specific 
material  system  need  be  established.  Hence  several  preliminary  MD  simulation  analysis  runs 
were  used  to  test  and  obtain  the  appropriate  dynamics  analysis  time,  appropriate  pressure  and 
temperature  controls  to  be  used  as  well  as  the  appropriate  material  system  size. 

3.4.1  Dynamic  time  average.  As  discussed  in  the  background  sections,  MD  simulations 
are  based  on  time-averaged  values  of  the  computed  thermodynamic  and  physical  properties  of 
the  molecular  system.  The  longer  the  time  dynamic  duration  of  a  MD  simulation,  the  better  the 
time  average  would  be.  However,  in  order  to  make  effective  use  of  the  available  and  limited 
computational  resources,  an  appropriate  duration  of  the  simulation  must  be  established. 

3.4.1. 1  Determination  of  dynamic  time.  The  dynamic  simulation  time  was  chosen  based 
on  the  dynamics  time  duration  of  100  ps,  200  ps  and  300  ps.  In  order  to  choose  an  appropriate 
dynamic  time  for  the  C-S-H  Jennite  material  system,  the  lxlxl  unit  cell  of  the  traditional  Jennite 
molecular  structure  was  subjected  to  energy  relaxation  using  COMPASS  forcefield  and  NVT 
ensemble  using  different  temperature  controls.  The  dynamics  time  was  obtained  based  on  the 
convergence  of  the  energy. 

Figure  12,  shows  the  energy  fluctuation  during  the  dynamics  analysis  for  these  different 
dynamic  times  durations  with  Nose  thermostat.  As  shown  in  figure  12,  there  was  a  good 
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convergence  in  the  energy  values  for  all  these  three  dynamic  durations.  Figures  12  and  13 
(simulations  run  using  Andersen  and  Velocity  Scale  thermostat)  also  showed  good  convergence 
for  all  the  three  transient  dynamic  analysis  durations.  However,  with  Berendsen  thermostat  as 
shown  in  figure  14,  a  longer  dynamic  time  was  required  for  the  energy  convergence.  Based  on 
the  above  simulation  results,  in  other  to  provide  a  good  accuracy  and  computational  efficiency, 
the  dynamic  duration  of  200  ps  was  chosen  for  all  the  full  analysis  studies.  In  the  case  of 
Magnesium  modified  Jennite,  the  systematic  study  of  the  dynamics  time  duration  was  performed 
with  Universal  Force  Field  and  was  found  to  show  similar  behavior  as  in  the  case  of  traditional 
C-S-H  Jennite.  See  Appendix  A1  for  the  NVT  dynamics  for  the  data  on  54.8%  Mg  modified 
Jennite. 
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Figure  12.  Total  energy  versus  simulation  time  -  Nose  thermostat 
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Figure  13.  Total  energy  versus  simulation  time  -  Andersen  thermostat 
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Figure  14. Total  energy  versus  simulation  time  -  Velocity  Scale  thermostat 
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Figure  1 5. Total  energy  versus  simulation  time  -  Berendsen  thermostat 

3.4.2  Thermodynamic  temperature  and  pressure  control.  The  thermodynamic 
parameter  of  temperature  and  pressure  in  MD  simulations  is  based  on  kinetic  energy  of  the  atom 
configurations.  A  stable  temperature  control  is  obtained  by  different  computational  temperature 
control  methods.  A  choice  needs  to  be  made  from  the  various  temperature  and  pressure  control 
algorithms  that  was  discussed  in  section  2.3.  In  order  to  do  this,  MD  system  equilibration  using 
NVT  and  NPT  ensembles  were  used. 

3.4.2. 1  Determination  of  temperature  control  method.  To  determine  the  temperature 
control  algorithm  to  be  used,  the  NVT  dynamic  runs  discussed  earlier  {figures  12  to  15)  were 
used.  From  these  figures,  the  dynamic  energy  fluctuations  in  Nose  thermostat  tends  to  be  smooth 
and  less  noisy  than  the  Andersen  and  Velocity  Scale  thermostat.  The  Berendsen  result  showed  a 
larger  fluctuation  and  is  noisy,  and  was  not  selected. 
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3.4.2.2  Determination  of  dynamic  pressure  control  method.  The  already  established 
dynamic  time  of  200ps  and  Nose  thermostat  were  used  to  determine  the  pressure  control  method 
to  use.  An  NPT  equilibration  dynamics  based  on  a  single  unit  cell  traditional  Jennite,  using 
COMPASS  forcefield  was  conducted  with  different  pressure  control  algorithms.  The  results  are 
shown  in  the  figures  16  -18. 
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Figure  16.  Total  energy  versus  simulation  time  -  Parinelo  barostat 


It  is  seen  that  using  the  Parinello  barostat  provides  a  smoother  and  less  noisy  total  energy 


profile  during  pressure  equilibration  than  those  from  the  Andersen  and  Berendsen  barostat. 
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Figure  1 7.  Total  energy  versus  simulation  time  -  Andersen  barostat 
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Figure  18.  Total  energy  versus  simulation  time  -  Berendsen  barostat 
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3.4.3  Material  system  size.  The  MD  simulation  uses  periodic  boundary  conditions  in 
order  to  minimize  the  boundary  effects  on  boundary  atoms,  and  simulate  a  larger  material 
configuration.  However,  a  critical  system  size  that  shows  a  convergence  of  the  physical 
properties  should  be  established  for  a  given  molecular  system.  In  the  present  study,  we  seek  a 
material  system  size  that  gives  a  good  convergence  on  elastic  stiffness  modulus  and  total 
energies  as  the  material  system  size  increases.  The  appropriate  cell  size  was  chosen  from  the 
following  possibilities  using  the  MD  simulations  conducted  for  these  sizes:  1  unit  cell  (lxlxl),  8 
unit  cells  (2x2x2),  27  unit  cells  (3x3x3),  64  unit  cells  (4x4x4),  and  125  unit  cells  (5x5x5). 

3.4.3. 1  Determination  of  material  system  size.  In  addition  to  a  single  unit  cell 
configuration,  8units,  27  units,  64  units  and  125  units  supercells  were  created  from  single  unit 
cell  traditional  Jennite  and  one  of  the  modified  Jennite  (54.8%Mg)  in  the  present  work.  Figure 
19  shows  the  pictures  of  1,  8  and  64  unit  cells  of  54.8%  modified  Jennite  structure. 


luc  8uc  64uc 

Figure  19.  Material  system  sizes,  1,  8  and  64  unit  cells  of  a  Mg  modified  Jennite 

All  these  different  configurations  for  both  traditional  C-S-H  Jennite  and  Magnesium 
modified  C-S-H  Jennite  were  employed  for  the  MD  analysis.  The  time  averaged  total  energies 


43 


and  the  elastic  moduli  obtained  were  compared.  Figures  20  and  21  show  that  there  is  a 
convergence  of  the  average  total  energies  per  Mol  of  the  structure  as  the  size  increases. 
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Figure  20.  Normalized  average  total  energy  for  traditional  C-S-H  structures:  1,  8,  27,  64,  125 
unit  cells 
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Figure  21.  Normalized  average  total  energy  for  55%Mg  modified  structures:  1,  8,  27,  64,  125 
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Likewise  the  elastic  stiffness  modulus,  determined  using  the  Hill  average,  shows  that  as 
the  size  of  the  cell  increases  the  moduli  converges  for  both  traditional  and  modified  Jennite.  See 
figures  22  and  23. 
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Figure  22.  Elastic  stiffness  modulus  (Hill  average)  for  traditional  C-S-H  structures:  1,  8,  27,  64 
and  125  unit  cells. 
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Figure  23.  Elastic  stiffness  modulus  (Hill  average)  for  55%Mg  Modified  C-S-H  structures:  1,  8, 


27  and  64  unit  cells 
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From  figures  20-23,  the  system  size  of  64  unit  cells  was  chosen  as  the  size  of  the  material 
molecular  structure  configuration  to  work  with  in  future  MD  analysis  simulations  for  the 
material  modeling.  This  material  system  size  showed  good  convergence  and  is  within  the 
computational  capacity  available  to  us  in  the  present  time.  This  molecular  system  size  will  be 
sufficient  to  handle  the  dynamics  analysis  within  reasonable  computational  time,  while  giving 
reliable  results  to  understand  the  effect  of  magnesium  ion  exchanges. 

3.4.4  Influence  of  forcefield.  As  discussed  earlier  on  the  background  of  MD,  the  choice 
of  forcefield  plays  an  important  role  in  the  fidelity  of  results  obtained  from  MD  simulations. 
Hence  in  order  to  choose  between  the  UNIVERSAL  and  COMPASS  forcefields,  a  NVT 
dynamic  simulation  was  performed  on  two  of  modified  and  traditional  Jennite  structures  using 
the  64  unit  cell  size,  Nose  thermostat  and  200ps  dynamics  analysis  time,  results  obtained  is 
shown  in  figure  24. 


Figure  24.  Total  energy  profile  for  64  unit  cell  traditional  and  modified  Jennite  (55%Mg)  using 


COMPASS  and  universal  forcefields. 
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From  figure  24,  the  values  of  the  total  energy  for  the  traditional  Jennite  and  the  55% 
modified  Jennite  when  Universal  forcefield  was  used  is  higher  than  the  values  obtained  when 
COMPASS  forcefield  was  used.  Other  modified  structures  were  also  studied  using  the  same  MD 
simulation  analysis  with  both  forcefields.  Similar  result  is  obtained  when  33%  modified  Jennite 
( figure  25)  and  83%  modified  Jennite  {figure  26)  were  used. 


xio5 


Figure  25.  Total  energy  profile  for  64  unit  cell  traditional  and  modified  Jennite  (33%Mg)  using 
COMPASS  and  UNIVERSAL  forcefields. 

In  all  the  cases,  it  was  also  observed  that  the  COMPASS  forcefield  profile  appears 
smoother  with  less  noise  fluctuations  when  both  forcefield  were  plotted  on  the  same  scales 
figures  24-26.  Based  on  these,  COMPASS  forcefield  was  chosen  and  used  in  the  material 
property  determinations  as  discussed  in  Chapter  4. 
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Figure  26.  Total  energy  profile  for  64  unit  cell  traditional  and  modified  Jennite  (83%Mg)  using 
COMPASS  and  universal  forcefields. 

3.4.5  MD  parameters  summary.  Following  the  above  discussions,  in  summary,  the 
parameters  used  for  the  mechanical  property  determination  and  shear  deformation  studies  to 
understand  the  influence  of  Magnesium  ion  exchanges  in  the  current  work  are; 

•  Dynamic  time:  200  ps 

•  Thermostat:  Nose 

•  Barostat:  Parinello 

•  Cell  size:  64  unit  cell 


Forcefield:  COMPASS 
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CHAPTER  4 


Results  and  Discussions 

The  focus  of  the  present  work  is  to  understand  the  effect  of  magnesium  ion  exchange  in 
traditionally  calcium  based  C-S-H  on  mechanical  stiffness  properties  and  the  shear  deformation 
behavior  following  the  material  modeling  methodology  outlined  in  chapter  2.  These  are  now 
presented  and  discussed  in  this  chapter. 

4.1  Elastic  Stiffness  Modulus 

The  MD  modeling  analysis  parameters  obtained  and  discussed  in  chapter  3  were  used  to 
investigate  the  influence  of  the  Magnesium  inclusion  and  ion  exchange  using  the  64  unit  cell 
structure  of  C-S-H  Jennite  structure  at  a  temperature  of  298K.  The  elastic  modulus,  mechanical 
properties  for  the  traditional  calcium  based  and  the  various  magnesium  based  molecular 
configurations  of  C-S-H  Jennite  were  obtained.  Figure  27  shows  the  initial  and  final  molecular 
structures  -  after  dynamics  simulations  on  the  systems. 


JO:  Traditional  Jennite  (0%Mg) 


Jl:  Modified  Jennite  (33%Mg) 


J3:  Modified  Jennite  (55%Mg) 


J4:  Modified  Jennite  (83%Mg) 


J5:  Modified  Jennite  (100%Mg) 


Figure  27.  Initial  and  final  Jennite  structures  -  Traditional  and  Mg  Ion  Exchanged  C-S-H 
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Figure  28  shows  the  average  total  energy  for  various  magnesium  modified  C-S-H  structures  as  a 
function  of  various  magnesium  percentages  in  the  modified  C-S-H  structures. 
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Figure  28.  Change  in  total  energy  per  mole  for  various  Mg  ion  exchange  percentages 
The  elastic  moduli  was  determined  using  Hill’s  averages  for  the  traditional  and 
Magnesium  ion  exchanged  64  unit  cell  C-S-H  Jennite  structures.  Figure  29  shows  the  variation 
in  the  predicted  elastic  moduli  with  COMPASS  forcefield  following  the  MD  based  material 
modeling  methodology. 


Latter  studies,  involving  Calcium  ion  replacement  of  all  1  to  9  calcium  atoms  was  also 
performed  to  obtain  mechanical  stiffness  properties  for  other  calcium  replacement  percentages  of 
7%,  23%,  43%  ,  68%Mg,  and  100%  Mg.  In  addition  to  this;  statistical  analysis  with  95% 
confidence  level  was  performed  to  obtain  the  number  of  different  variations  of  each  modified 
structures  (with  respect  to  location  of  the  calcium  ion  replaced)  that  is  required  to  have  a  95% 
confidence  level,  and  was  analyzed.  Appendix  A2  shows  the  table  of  all  variants  obtained. 
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Figure  29.  Predicted  elastic  modulus  (Hill  Average)  for  various  Mg  ion  exchange  percentages  of 
64  unit  cell  C-S-H  Jennite  structures-COMPASS  forcefield. 

Figure  30  shows  the  mean  values  of  the  normalized  average  total  energies  for  all  the  Mg 
ion  exchange  percentages  in  C-S-H  Jennite  structure.  The  errors  associated  with  the  energies  are 
very  small  and  not  visible  form  figure  30.  Appendix  A3  show  the  table  of  the  result  of  the  total 
energies  and  the  corresponding  error  associated  with  each  Mg  modified  structure  following  the 
statistical  analysis  of  the  atom  locations  employed  for  the  magnesium  exchange. 

In  the  same  manner,  the  various  variants  of  the  Mg  modified  Jennite  were  also  used  to 
obtain  the  elastic  stiffness  modulus,  and  figure  31  shows  the  mean  values  of  the  elastic  modulus 
versus  the  %Mg. 
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Figure  30.  Total  energy  per  mole  for  various  Mg  ion  exchange  percentages. 


Figure  31.  Predicted  elastic  modulus  (Hill  Average)  for  various  Mg  ion  exchange  with  the 
variants  and  associated  errors-COMPASS  forcefield. 

The  variation  in  the  predicted  elastic  modulus  values  presented  in  figure  29  shows  a  sharp 
increase  in  the  elastic  modulus  of  Jennite  from  the  traditional  for  33%  magnesium  ion  exchange, 
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followed  by  reduction  in  the  elastic  modulus  as  the  percentage  of  the  magnesium  increases 
further. 

The  values  of  the  predicted  mechanical  and  other  material  properties  depend  on  the 
forcefield  that  defines  the  energy  of  the  molecular  level  interactions.  A  slightly  different  result 
with  higher  values  was  observed  when  universal  forcefield  was  used  to  determine  the  elastic 
moduli.  Figure  32  shows  the  variation  in  the  predicted  modified  elastic  modulus  for  different 
magnesium  ion  exchange  percentages  using  universal  forcefield. 
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Figure  32.  Predicted  elastic  modulus  (Hill  Average)  for  various  Mg  ion  exchange  percentages  of 
64  unit  cell  C-S-H  Jennite  structures-Universal  forcefield. 

The  corresponding  predicted  bulk  and  shear  modulus  and  the  as  well  as  the  Poisson’s 
ration  were  obtained  and  are  presented  in  figures  33-35  below.  Negative  values  obtained  for  the 
predicted  bulk  moduli  and  Poisson’s  ratios  are  based  on  the  mechanics  based  theories  of  defined 
effective  bulk  modulus  from  compliance  matrix  as  discussed  in  Chapter  2  and  needs  further 
investigation. 
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Figure  33.  Predicted  shear  modulus  (Hill  Average)  variations  for  various  magnesium  ion 
exchange  percentages  studied. 
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Figure  34.  Predicted  bulk  modulus  (Hill  Average)  variations  for  various  magnesium  ion 
exchange  percentages  studied. 
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Figure  35.  Predicted  Poisson  ratio  (Hill  Average)  variations  for  various  Magnesium  ion 
exchange  percentage  studied. 

4.1.1  Mechanical  properties  summary.  Elastic  modulus  investigation  showed  that 
material  chemistry  does  affect  structural  material  stiffness  modulus  properties,  and  varies  as  the 
amount  of  Mg  inclusion  increases.  These  predicted  material  properties  are  based  on  MD 
simulations  and  material  chemistry  level  modeling  and  thus  provides  an  effective  computational 
methodology  to  understand  the  expected  variations  due  to  material  chemistry  changes  such  as 
due  to  Magnesium  ion  exchange  of  C-S-H  Jennite.  Further  investigation  is  needed  to  explain 
some  of  the  non-physical  computational  modeling  values  noted  in  the  predicted  values  of  bulk 
modulus  and  Poisson’s  ratio  that  are  based  on  mechanics  based  approximations  following  the 
compliance  matrix  definitions  in  MD  simulation  analysis. 

The  prior  MD  material  modeling  predictions  are  based  on  the  time  averaged  values  of  the 
modulus  values  from  derivatives  of  energy  as  defined  in  chapter  2  for  effective  values  of 
material  stiffness  modulus.  Engineering  scale  properties  of  interest  is  the  stress-strain 
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deformation  behavior  and  effective  prediction  of  estimated  stress-strain  behavior  due  to  material 
chemistry  changes  in  the  systems,  such  as  those  discussed  due  to  Magnesium  exchange  in  C-S-H 
Jennite.  In  the  present  work,  the  influence  of  magnesium  exchange  on  shear  deformation  of  C-S- 
H  Jennite  molecular  structure  is  studied.  Details  of  the  shear  deformation  with  emphasis  on  the 
Magnesium  ion  exchanged  C-S-H  Jennite  is  presented  next. 

4.2  Shear  Deformation 

Shear  deformation  was  implemented  by  imposing  a  constant  volume  angular  deformation 
change  on  each  crystallographic  plane  of  the  triclinic  C-S-H  Jennite  structure  and  the 
Magnesium  modified  systems.  In  the  shear  deformation  process,  the  C-S-H  molecular  structures 
were  subjected  first  to  NPT  conditions.  Then,  they  were  subjected  to  an  incremental  angular 
change  making  sure  the  volume  remains  constant  defining  the  shear  deformation  process  as 
outlined  in  chapter  2.  The  resulting  structure  was  relaxed  after  which  an  NVT  dynamic  analysis 
for  each  of  the  deformed  structural  configuration  was  performed  to  obtain  the  corresponding 
stress  value  following  the  Virial  definition  of  shear  stress.  The  shear  deformation  process  was 
applied  to  all  three  triclinic  faces  of  C-S-H  Jennite  as  shown  in  figure  36. 
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Figure  36.  Shear  deformation  modeling  along  three  crystallographic  planes  of  triclinic  C-S-H 
structure. 

The  schematics  presented  Figure  3  7  shows  the  procedure  of  the  MD  modeling  analysis 
for  the  shear  deformation.  To  define  the  complete  shear  deformation  process,  the  original  and 
quasi-static  deformation  of  the  triclinic  structure  was  performed  for  75  different  sheared 
configurations.  Each  shear  deformed  configuration  provides  the  corresponding  shear  stress, 
associated  strain  for  each  of  the  deformed  configuration  due  to  shearing  along  the  three  surfaces 


of  the  triclinic  traditional  and  modified  C-S-H  structures. 
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Figure  3  7.  MD  implementation  of  shear  deformation.  Showing  an  NPT  relaxation  of  the  initial 
structure,  followed  by  an  incremental  angular  change  -  the  deformation  step,  corresponding  NVT 
equilibration  of  the  deformed  structure  to  obtained  the  shear  stress  of  the  deformed  structure  and 
continuation  of  the  iterative  process. 

In  the  present  work,  two  of  the  Magnesium  ion  exchange  modified  structures  were 
analyzed  (33%Mg  and  83%Mg)  in  addition  to  the  traditional  Jennite  which  had  been  obtained 
earlier  in  our  research  group  [50].  Shear  stress-strain  deformation  on  the  three  shearing  planes 
cb,  ba  and  ac  of  the  triclinic  C-S-H  system  were  studied  and  compared. 
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4.2.1  Shear  deformation,  cZ>-plane.  Figures  38-40  shows  the  MD  modeling  predicted 
shear  stress-strain  curves  for  the  shear  deformation  along  the  cb-plane  for  the  traditional  C-S-H, 
33%  magnesium  modified  C-S-H  and  83%  magnesium  modified  C-S-H  Jennite  structures 
respectively. 


Traditional  C-S-H,  (0%Mg) 


Figure  38.  Predicted  Shear  Stress-strain  deformation  behavior  for  shearing  along  the  cb-plane  for 
traditional  Jennite.  Ref  [50] 

Employing  the  predicted  stress-strain  deformation  behavior,  shear  modulus  was  obtained, 
as  the  slope  of  the  linear  region  of  the  shear  stress  -  strain  deformation  behavior  shown  in  figure 
38-40.  The  ultimate  shear  stress  was  defined  as  the  maximum  stress  observed  within  the  shear 


strain  of  the  shearing  deformation  studied. 
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Modified  C-S-H,  (33%Mg) 


Figure  39 .  Predicted  Shear  Stress-strain  deformation  behavior  for  shearing  along  the  c/i-plane 
for  33%  magnesium  modified  Jennite 


Modified  C-S-H,  (83%Mg) 


Figure  40.  Predicted  Shear  Stress-strain  deformation  behavior  for  shearing  along  the  r7?-plane  for 
83%  magnesium  modified  Jennite 
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The  MD  material  modeling  predicted  values  for  shear  modulus,  strength  and  ultimate 
shear  stress  is  presented  in  table  5.  These  engineering  material  properties  for  traditional  C-S-H 
and  magnesium  modified  Jennite  have  been  solely  obtained  from  MD  modeling  analysis  of  the 
material  chemistry  structure  and  provides  an  effective  material  modeling  methodology  to 
understand  the  estimated  variations  in  the  stress  -  strain  deformation  behavior  due  to  material 
chemistry  changes. 

Table  5 

Shear  Modulus,  Strength  and  Maximum  Shear  Deformation  of  cb-plane 


Plane-  CB 

0%Mg 

33%Mg 

83%Mg 

Shear  Modulus  (GPa) 

1 1.0±2.4 

12.88±  1.21 

18.19  ±  1.75 

Shear  Strength  (GPa) 

0.5193 

0.8014 

1.038 

Ultimate  Stress  (GPa) 

1.1 

1.275 

1.48 

The  plot  of  the  shear  moduli  and  the  ultimate  shear  stress  is  also  presented  in  figures  41 
and  42.  It  is  observed  that  as  the  percentage  of  Magnesium  increases,  these  two  predicted 
material  properties  also  show  an  increase  for  shear  deformation  of  c6-plane. 
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Shear  deformation  of  the  traditional  and  Magnesium  modified  C-S-H  Jennite  along  other 
crystallographic  planes  is  presented  next. 

4.2.2  Shear  deformation,  Z>a-plane.  Shear  stress  -  strain  deformation  behavior  due  to 
shearing  along  the  crystallographic  6a-plane  is  presented  in  figures  43-45  for  traditional  and 
Magnesium  modified  C-S-H. 


Traditional  C-S-H 


Figure  43.  Predicted  Shear  Stress-strain  deformation  behavior  for  shearing  along  the  ba- plane 


for  traditional  Jennite.  Ref  [50] 
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Modified  C-S-H,  (33%Mg) 


Figure  44.  Predicted  Shear  Stress-strain  deformation  behavior  for  shearing  along  the  ha-plane 
for  33%  magnesium  modified  Jennite 


Modified  C-S-H,  (83%Mg) 


Figure  45.  Predicted  Shear  Stress-strain  deformation  behavior  for  shearing  along  the  ha-plane 
for  83%  magnesium  modified  Jennite 
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The  corresponding  shear  moduli  and  ultimate  shear  stress  obtained  for  shear  deformation  along 
/v/ -plane  is  presented  in  figures  46  and  47. 


Figure  46.  Shear  Modulus  variations  due  to  Mg  Ion  Exchange,  6n-plane 
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Figure  47.  Ultimate  shear  stress  variations  due  to  Mg  Ion  Exchange,  /w-plane 
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4.2.3  Shear  deformation,  ac-plane.  Shear  stress  -  strain  deformation  behavior  due  to 
shearing  along  the  crystallographic  ac-plane  is  presented  in  figures  48-50  for  traditional  and 
Magnesium  modified  C-S-H. 


Traditional  C-S-H,  (0%Mg) 


Figure  48.  Predicted  Shear  Stress-strain  deformation  behavior  for  shearing  along  the  ac-plane  for 


traditional  Jennite.  Ref  [50] 
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Modified  C-S-H,  (33%Mg) 


Figure  49.  Predicted  Shear  Stress-strain  deformation  behavior  for  shearing  along  the  ac-plane  for 
33%  magnesium  modified  Jennite. 


Modified  C-S-H,  (83%Mg) 


Figure  50.  Predicted  Shear  Stress-strain  deformation  behavior  for  shearing  along  the  ac-plane  for 
83%  magnesium  modified  Jennite 
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The  corresponding  shear  moduli  and  ultimate  shear  stress  obtained  for  shear  deformation  stress  - 
strain  curves  along  ac-plane  is  presented  in  figures  51  and  52. 


Figure  51.  Shear  Modulus  variations  due  to  Mg  Ion  Exchange,  ac-plane 


Figure  52.  Ultimate  shear  stress  variations  due  to  Mg  Ion  Exchange,  ac-plane 
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4.2.4  Shear  deformation  analysis  results  summary.  The  results  discussed  earlier  from  MD 
material  modeling  methodology  for  shear  deformation  along  the  three  crystallographic  planes  of 
triclinic  C-S-H  structure  show  the  following  changes  in  the  predicted  material  properties. 

•  Predicted  shear  modulus  for  shear  deformation  along  all  three  crystallographic  planes 
increases  as  the  %Mg  increases 

•  The  predicted  ultimate  shear  strength  before  material  failure  increases  as  the  %Mg 
increases 

Table  6  summarizes  the  predicted  shear  modulus  values  for  shear  deformation  along  the  three 
crystallographic  planes  of  traditional  and  Magnesium  ion  exchange  modified  C-S-H 
configuration  studied.  The  corresponding  ultimate  shear  strength  obtained  from  the  shear 
deformation  behavior  along  the  three  crystallographic  planes  of  triclinic  C-S-H  molecular  system 
is  presented  in  Table  7. 

Table  6 


Predicted  Shear  Modulus  (GPa)  all  planes 


Shear  Modulus  (GPa) 

%  Mg 

0% 

33% 

83% 

ac-Plane 

9.7±1.1 

15.66±1.9 

16.91±0.7 

cb-Plane 

11.0±2.4 

12.88±1.2 

18.19±1.75 

/v/- Plane 

13.1±1.0 

18.38±0.7 

27.55±3.6 
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Table  7 


Ultimate  Shear  Stress  (GPa)  all  planes 


Ultimate  Shear  Stress  (GPa) 

%  Mg 

0% 

33% 

83% 

Plane  -  ac 

1.2 

1.302 

1.287 

Plane  -  cb 

1.1 

1.275 

1.48 

Plane  -  ba 

1.3 

1.536 

2.038 

Next  section  discusses  a  closer  look  at  the  shear  deformation  process  and  the  atom  by 
atom  deformation  study.  The  traditional  Jennite  and  two  modified  Jennite  structure  were 
analyzed  below  and  is  required  to  understand  the  shear  deformation  and  failure  captured  in  the 
stress-strain  curves. 

4.3  Shear  Deformation  Atom  Trajectory  Study 

The  shear  deformation  response  was  studied  by  tracking  the  positions  of  the  atoms  as  the 
deformation  proceeds  and  to  observe  the  movement  of  each  atom  as  the  stress  is  applied  to  the 
structure.  In  order  to  do  this,  the  atoms  in  the  initial  structure  were  grouped  into  10  different 
sections  or  layers.  C-S-H  molecular  structure  shows  calcium  oxide  layers  with  free  water  and 
silica  chains  interspersed  between  these  layers.  Each  group/section  is  designated  as  layers  1  to 
10,  with  the  topmost  layer  being  the  tenth  layer.  The  location  of  CaO  molecules  was  used  to 
obtain  the  centroid  of  the  associated  CaO  (Ca)  molecules.  The  locations  of  the  centroids  of  the 
CaO  atoms  in  a  single  layer  were  consequently  used  to  determine  the  centroid  of  each  layer.  In 
the  present  analysis,  the  centroid  is  independent  of  the  atomic  mass  of  the  constituent  atoms.  So 
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that  when  Ca  is  replaced  with  Mg,  the  effect  of  mass  is  not  taken  into  account.  The  displacement 
each  centroid  is  tracked,  as  the  material  is  deformed.  The  actual  displacement  of  the  centroid 
versus  strain  curves  was  obtained  for  each  layer.  Also  a  hypothetical  linear  elastic  behavior 
curve,  which  assumes  that  the  material  continues  to  behave  elastically  throughout  the  strain 
region  was  also  plotted  for  each  layers.  This  was  done  for  the  traditional  Jennite  and  the  33%  and 
83%  magnesium  modified  Jennite  for  all  the  crystallographic  planes.  The  strain  region  on  the 
displacement-strain  curves  where  there  is  non-linear  movement  of  the  layer  centroid  was 
correlated  to  the  strain  region  when  non-linear  deformation  starts  to  manifest  in  their  respective 
stress-strain  curves.  The  results  of  6<:/ -planes  is  compared  and  presented  next. 

4.3.1  Shear  deformation,  atom  trajectory  study,  ba -plane.  The  movement  of  the 
centroid  of  each  layer  along  the  x-axis  as  the  strain  is  increased  is  presented.  Figures  53  and  54 
show  the  displacement-strain  curves  along  the  ba- plane,  for  layers  3-6  and  layers  7-10 
respectively  for  the  traditional  Jennite.  Figures  56  and  57  shows  the  displacement-strain  curves 
for  33  %  modified  Jennite  along  the  ba-plane.  Also,  Figures  59  and  60  show  those  for  the  83% 
modified  Jennite,  along  the  6c/-plane.  The  figures  (53  and  54)  also  show  the  hypothetical 
continued  linear  displacement  behavior  of  each  material  with  strain  as  they  undergo  shear 
deformation.  The  hypothetical  curve  assumes  that  the  material  behaves  linearly  throughout  the 
strain  region.  The  actual/observed  behavior  is  different  from  the  hypothetical  behavior,  the 
deformations  in  most  cases  follow  linear  deformation  but  as  the  strain  level  increases,  they  began 
to  exhibit  deviation  from  the  linear  behavior.  The  strain  values  when  there  is  a  significant  change 
in  linear  to  non-linear  behavior,  depicted  by  either  the  change  in  slope  or  deviation  from  the 
hypothetical  curve  is  observed  and  compared  with  the  strain  values  at  highest  stress  in  the  stress- 


strain  curves. 
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Traditional  Jennite,  ba-Plane 


(c)  Layer-5  Strain  (d)  Layer-6  Strain 


Figure  53.  Atomic  Centroid  Displacement-Strain  Curves  for  Layers  3-6  along  x-axis  for 
Traditional  Jennite  when  sheared  along  ba- plane. 

From  Figures  53  and  54,  and  as  expected,  the  first  few  bottom  layers  (3-5)  appear  to 
move  only  slightly  along  the  x-axis  compared  to  the  upper  layers  (7-9)  for  all  the  Jennite  types. 
The  atomic  centroid  movement  of  the  first  four  layers  for  all  the  Jennite  appears  to  move  slightly 
less  than  the  centroids  of  the  topmost  six  layers. 
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Traditional  Jennite,  ba-Plane 


(c)  Layer-9  Strain  (d)  Layer-10  Strain 


Figure  54.  Atomic  Centroid  Displacement-Strain  Curves  for  Layers  7-10  along  x-axis  for 
Traditional  Jennite  when  sheared  along  ha-plane 

The  expanded  inset  in  the  figures  shows  vertical  lines  drawn  through  the  location  where 
the  non-linear  behavior  is  found.  This  is  found  to  be  strain  of  about  0.04  for  layer-3.  The  strain 
values  for  layers  4-6  in  Figure  53  are  0.05,  0.055,  and  0.075,  while  those  for  layers  7-10  in 
Figure  54  are:  0.065,  0.068,  0.077  and  0.07.  From  the  stress-strain  curve  of figure  55,  the 
traditional  Jennite  experiences  the  highest  stress  and  transition  to  non-linear  behavior  within  the 
strain  range  of  0.04  to  about  0.076,  this  values  corresponds  to  the  strain  region  in  the 
displacement-strain  curves  {Figures  53  and  54)  where  non-linear  displacement  of  the  atom  is 


found. 
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Figure  55.  Stress-Strain  Curve  for  Traditional  Jennite  when  sheared  along  the  /?<:/- Plane 

This  displacement  of  the  atom  hence  gives  the  region  where  the  material  starts  to  yield 
and  shows  a  linear  to  non-linear  transition  due  to  the  applied  deformation. 

Similarly  for  the  33%Mg  modified  Jennite,  the  non-linear  deformation  starts  occurring  at 
a  strain  between  (0.06,  0.055,  0.048,  0.074,  0.07,  0.047,  0.056  and  0.068)  for  layers  3-10 
respectively  as  depicted  in  Figures  56  and  57  below. 
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33%  Mg  Modified  Jennite,  ba-Plane 


(c)  Layer-5  strain  (d)  Layer-6  strain 


Figure  56.  Atomic  Centroid  Displacement-Strain  Curves  for  Layers  3-6  along  x-axis  for  33% 
Magnesium  Modified  Jennite  when  sheared  along  ba- plane 


The  strain  values  from  the  displacement-strain  curves  above  were  also  found  to 
correspond  to  the  strain  at  the  yield  region  in  the  stress-strain  curve  (0.04-0.076)  for  the  /w-plane 
of  the  33%  modified  Jennite  as  shown  in  Figure  58. 
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33%  Mg  Modified  Jennite,  ba-Plane 


(c)  Layer-9  Strain  (d)  Layer-10  Strain 


Figure  57.  Atomic  Centroid  Displacement-Strain  Curves  for  Layers  7-10  along  x-axis  for  33% 
Magnesium  Modified  Jennite  when  sheared  along  Z>«-plane 
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Figure  58.  Stress-Strain  Curve  for  33%  Mg  Modified  Jennite  when  sheared  along  the  6a- Plane 
Lastly,  83%  traditional  Jennite  also  experience  a  similar  strain  values,  first  there  was 
slight  non-linear  deformation  around  values  of  0.035-0.05  in  most  of  the  layers.  However  there 
was  a  return  to  linearity  before  a  final  non-linear  deformation  at  a  later  strain  values  of;  0.055, 
0.06,  0.057,  0.07,  0.064,  0.068,  0.07  and  0.08  for  layers  3  -10  respectively  ( Figures  59  and  60). 
The  same  observation  can  be  noticed  in  the  stress-strain  curve  as  depicted  in  Figure  61. 
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Figure  59.  Atomic  Centroid  Displacement-Strain  Curves  for  Layers  3-6  along  x-axis  for  83% 
Magnesium  Modified  Jennite  when  sheared  along  ha-plane 
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83%  Mg  Modified  Jennite 


Figure  60.  Atomic  Centroid  Displacement-Strain  Curves  for  Layers  7-10  along  x-axis  for  83% 
Magnesium  Modified  Jennite  when  sheared  along  ba- plane 
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Figure  61.  Stress-Strain  Curve  for  83%  Mg  Modified  Jennite  when  sheared  along  the  /w-Plane 
In  general,  the  strain  values  in  the  stress  strain  curve  falls  between  0.04-0.076,  which  is 
similar  to  those  from  the  displacement  strain  curve  for  all  the  layers  considered.  From  the 
observations,  it  can  be  seen  that  the  movement  of  the  atoms  can  give  information  about  the  strain 
region  where  the  material  starts  behaving  in-elastically  and  a  transition  from  linear  to  non-linear 
behavior  occurs. 


The  same  analyses  were  done  for  cb- plane  and  ac-plane,  the  displacement-strain  curves 
for  layers  3-10  for  33%  Mg  modified  Jennite  for  c'6-plane  and  83%  Mg  modified  Jennite  for  ac- 
plane  can  be  found  in  Appendices  B 1  and  B2  respectively. 
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CHAPTER  5 

Discussion  and  Future  Research 

5.1  Discussion  of  Results 

From  the  results,  we  have  been  able  to  establish  that  material  chemistry  of  traditional  C- 
S-H  Jennite  will  be  affected  by  ionic  replacement  of  the  Ca  ions  by  Mg  ions.  By  including 
Magnesium  in  increasing  proportion  in  the  C-S-H  Jennite,  the  total  energy  of  the  structure  was 
observed  to  decrease  with  the  increase  in  the  amount  of  Magnesium  added.  By  carrying  out  the 
computational  mechanical  stiffness  property  evaluation  on  the  structures,  the  elastic  modulus  of 
the  system  structures  were  determined.  It  was  found  to  decrease  as  the  amount  of  Magnesium  in 
the  molecular  structure  increased.  COMPASS  force  field  and  UNIVERSAL  force  field  tends  to 
describe  the  same  behavior  of  the  total  energy  of  the  molecular  system,  although  the  values 
obtained  using  UNIVERSAL  force  field  were  found  to  be  higher.  It  should  be  noted  that 
experimental  determination  of  these  elastic  properties  is  a  complex  task.  What  we  have  been  able 
to  establish  here  is  that  material  chemistry  play  important  role  in  the  elastic  stiffness  properties 
of  the  C-S-H  Jennite  structure.  The  trends  have  been  established  by  using  the  two  force  fields, 
but  there  is  a  need  to  further  improve  the  fidelity  of  the  elastic  stiffness  modulus  with  further 
studies  and  comparative  verification  and  validations. 

Results  obtained  from  the  study  of  the  influence  of  Magnesium  on  shear  deformation, 
shows  that  by  deforming  the  traditional  and  modified  C-S-H  Jennite,  the  shear  modulus  obtained 
from  the  stress-strain  curves  shows  an  increasing  trend  as  the  amount  of  Magnesium  in  the 
system  increases.  We  refer  to  the  peak  stress  value  of  the  stress-stress  curves  as  the  ultimate 
shear  stress  that  can  be  absorbed  by  the  materials  within  0.025  strain,  and  any  further 
deformation  results  in  failure.  The  ultimate  shear  stress  values  were  seen  to  increase  as  the 
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amount  of  Magnesium  in  the  systems  increases.  We  can  also  establish  here  that  material 
chemistry  definitely  affects  the  engineering  scale  properties  of  the  cement  structure.  The  present 
results  are  based  on  an  assuming  that  structure  of  the  modified  C-S-H  Jennite  remains  triclinic  , 
even  though  Magnesium  ion  are  smaller  than  Calcium  ions.  Present  study  clearly  demonstrated 
that  computational  material  modeling  following  molecular  dynamics  analysis  methodology  is  an 
effective  way  to  predict  and  understand  the  effect  of  material  chemistry,  and  additive  changes  on 
the  stiffness  and  deformation  characteristics  in  cementitious  materials.  Such  methodology  is  also 
extendable  to  others. 

5.2  Recommendations  for  Future  Work 

Future  research  should  validate  the  results  obtained  in  the  present  work  by  using  the  other 
mineral  forms  of  C-S-H  Jennite  such  as  Tobermorite  14A°.  It  is  also  recommended  that  if 
computational  capacities  allows,  a  larger  system  size  should  be  used  to  validate  results.  Since  all 
results  were  obtained  with  material  chemistry  structure  and  computational  molecular  dynamics, 
future  attempts  should  investigate  potential  experimental  scale  relevant  validations  using  nano¬ 
indentation  and  investigate  potential  methods  for  low  length  scale  characterization  for  further 
validations  and  comparisons. 
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Appendix  A 1 

Energy  versus  Time  for  NVT  dynamics  of  1  unit  cell  modified  Jennite  (54.8%Mg)  using 
UNIVERSAL  forcefield. 
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Appendix  A2 

Spatial  location  of  Calcium  atom  in  Jennite. 
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The  number  of  sampled  variation  based  on  95%  confidence  level 
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Appendix  A3 


The  Magnesium  modified  Jennite,  the  variants  and  the  mean  averages  of  Total  Energies  and 
Elastic  stiffness  modulus  are  presented. 


J-14.8%Mg  --  64unit  cell  structure 

Variant  Name 

Ca  replaced 

Average  TE  (Kcal/Mol) 

E(GPa) 

J 1 5-i 

1,5 

-4212.42 

6.504627 

J15-j 

2,6 

-4211.7 

7.616457 

J15-k 

3,7 

-4207.3 

6.445543 

J15-1 

4,8 

-4209.48 

6.984129 

J15-m 

1,9 

-4211.5 

5.609928 

AVERAGE  =  -4210.48  6.632137 

STD  DEV  =  1.867129  0.66163 

J-(23.27%Mg)  --  64  unit  cell  structure 

Variant  Name 

Ca  replaced 

Average  TE  (Kcal/Mol) 

J23-i 

1,2,4 

-4250.69 

J23-j 

1,2,3 

-4247.81 

J23-k 

4,6,9 

-4251.55 

J23-1 

3,8,9 

-4243.57 

J23-m 

4,6,7 

-4241.46 

J23-n 

1,5,8 

-4245.6 

J23 -previous 

2,3,7 

-4246.56 

AVERAGE  =  -4246.75 

STD  DEV  =  3.363506 
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J-(32.7%Mg)  —  64  unit  cell  structure 

Variant  Name 

Ca  replaced 

Average  TE  (Kcal/Mol) 

E(GPa) 

J32-i 

1,2, 8, 6 

-4284.06 

40.43001 

J32-j 

3, 5, 7, 9 

-4274.06 

45.98485 

J32-k 

1,2, 4, 9 

-4278.11 

46.72645 

J32-1 

3, 5, 4, 8 

-4274.56 

38.22898 

J32-m 

6, 7, 1,5 

-4280.47 

41.11704 

AVERAGE  =  -4278.25  42.49747 

STD  DEV  =  3.737928  3.299879 

J-(43.21%Mg)  --  64  unit  cell  structure 

Variant  Name 

Ca  replaced 

Average  TE  (Kcal/Mol) 

E(GPa) 

J43-i 

1,2, 3, 4, 5 

-4304.306205 

J43-j 

1,3, 6, 8, 9 

-4292.287899 

J43-k 

2, 6, 7, 8, 9 

-4301.330568 

J43-1 

1,4, 5, 6, 7 

-4300.503784 

J43-m 

2, 3, 5, 6, 8 

-4303.081435 

J43-n 

1,3, 4, 7, 9 

-4302.097487 

19.60885334 

J43 -Previous 

2, 4, 5, 6, 8 

-4304.86911 

34.53704037 

AVERAGE  =  -4301.210927  27.07294686 

STD  DEV  =  3.916679392  7.464093511 
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J-(54.8%Mg)  —  64  unit  cell  structure 

Variant  Name 

Ca  replaced 

Average  TE(Kcal/Mol) 

E(GPa) 

J55-i 

1,2, 4, 7, 8, 9 

-4328.670814 

18.21657314 

J55-j 

1,2, 4, 6, 8, 9 

-4327.079064 

19.86393489 

J55-k 

2, 3, 5, 6, 7, 9 

-4318.831345 

19.43865362 

J55-1 

1,2, 6, 7, 8, 9 

-4329.254849 

13.69626296 

J55-m 

1,3, 4, 5, 7, 8 

-4327.839037 

17.24501927 

J55-n 

1,2, 4, 6, 7, 9 

-4328.162613 

24.59421684 

AVERAGE  =  -4326.63962  18.84244345 

STD  DEV  =  3.556345346  3.263649791 

J-(67.97%Mg)  --  64  unit  cell  structure 

Variant  Name 

Ca  replaced 

Average  TE  (Kcal/Mol) 

E(GPa) 

J68-i 

2, 3, 4, 6, 7, 8, 9 

-4348.358954 

15.55705204 

J68-j 

1,3, 4, 5, 7, 8, 9 

-4343.537117 

15.24017854 

J68-k 

1,2, 3, 4, 5, 6, 8 

-4348.308332 

12.6434633 

J68-1 

1,2, 4, 5, 7, 8, 9 

.4347.734238 

15.36978749 

J68-m 

1,2, 3, 5, 6, 7, 9 

-4347.420075 

12.65037346 

J6  8 -Previous 

1,2, 4, 6, 7, 8, 9 

-4345.533912 

18.05476625 

AVERAGE  =  -4346.815438  14.91927018 

STD  DEV  =  1.742156381  1.865084387 
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J-(82.9%Mg)  —  64  unit  cell  structure 

Variant  Name 

Ca  replaced 

Average  TE  (Kcal/Mol) 

E 

J83-i 

2,3,4,5,6,7,8,9 

-4361.005939 

10.2259845 

J83-j 

13,4,5,6,7,8,9 

-4361.090475 

14.53332378 

J83-k 

1,23,5,6,7,8,9 

-4359.934655 

14.75545451 

J83-1 

1,23,4,6,7,8,9 

-4360.238779 

15.53425973 

J83-m 

1,23,4,5,6,8,9 

-4359.450767 

14.1900116 

J8  3 -previous 

1,2,4,5,6,7,8,9 

-4361.676025 

14.3991761 

AVERAGE  =  -4360.47814  13.9397017 

STD  DEV  =  0.802252931  1.713958498 
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Appendix  B1 

Figures  showing  the  atomic  centroid  displacement-strain  curves  for  layers  3-10  along  for  33% 
Magnesium  Modified  Jennite  when  sheared  along  cb- plane. 
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Appendix  B2 

Figures  showing  the  atomic  centroid  displacement-strain  curves  for  layers  3-10  along  for  83% 
Magnesium  modified  Jennite  when  sheared  along  ac-plane. 
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Atomic  Centroid  Displacement-Strain  Curves  for  Layers  7-10  for  83%  Mg  modified  when 
sheared  along  ac-plane. 
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Appendix  Cl 

Shear  deformation  code 

#BIOSYM  btcl  3 

# 

#  Input  File  For  Discover  Generated  By  Materials  Studio 

#  Input  Client  Model  Document: 

#  Job:  Disco  Dynamics  # 
autoEcho  on 

# 

#  Begin  Forcefield  Section 
begin  forcefield  =  compass 

#  Nonbond  section: 

forcefield  nonbond  +separate_coulomb  \ 
vdw\ 

summation_method  =  ewald  \ 
ewald_accuracy  =  l.e-003  \ 
update_width  =  3.00  \ 
coulomb  \ 

summation_method  =  ewald  \ 
ewald_accuracy  =  l.e-003  \ 
update_width  =  3.00  \ 
dielectric_value  =  1.0000 

#  End  Forcefield  Section 


#Stage  Name:  minimize 
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minimize  \ 

iteration_limit  =  1000  \ 
execute  frequency  =  10  \ 
command  =  (print  output  +energy_summary) 

#  Dynamics  Section:  NPT  equilibration 

minimize 

dynamics  \ 

time  =  100000.00  \ 

timestep  =  1.00  \ 

ensemble  =  NPT  \ 

temperature  =  298.00  \ 

press_choice  =  stress  \ 
sxx  =  -0.00003  \ 
syy  =  -0.00003  \ 
szz  =  -0.00003  \ 
deviation  =  5000.00  \ 
execute  frequency  =  100.00  \ 

command  =  (print  table  file  =  uc64-SbaNPT.tbl  +average  +batch_average  batch_size  =100  +cell 
+energy  +state  +stress  +strain)  \ 
execute  frequency  =  1000.00  \ 

command  =  (print  archive  filename  =  uc64-SbaNPT.arc  +coordinates)  \ 
execute  frequency  =  1000.00  \ 

command  =  (writeFile  dynamics_restart  filename  =  uc64-SbaNPT.xdyn) 


writeFile  coordinate  filename  =  uc64-SbaNPT.car 
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#  Cell  parameters 

set  cP  list  [cellParameter] 
set  cP_a  [lindex  $cP_list  0] 

set  cP_b  [lindex  $cP _ list  1] 

set  cP_c  [lindex  $cP_list  2] 
set  cP_alpha  [lindex  $cP_list  3] 
set  cP  beta  [lindex  $cP_list  4] 
set  cP_gamma  [lindex  $cP_list  5] 

$a  =  $cP_a 
$b  =  $cP_b 
$c  =  $cP_c 
$alpha  =  $cP_alpha 
$beta  =  $cP_beta 
$gamma  =  $cP_gamma 

$pi  =  16.0*atan(l. 0/5.0)  -4*atan(l. 0/239.0) 

$Calp  =  cos($alpha*$pi/180) 

$Cbet  =  cos($beta*$pi/180) 

$Cgam  =  cos($gamma*$pi/180) 

$vol0  =  $a*$b*$c*sqrt(l-$Calp*$Calp-$Cbet*$Cbet-$Cgam*$Cgam+2*$Calp*$Cbet*$Cgam) 

#  Deformation  of  the  system 
for  (set  i  1)  ($i  <=  75)  (incr  i  1)  ( 

$gamma  =  $gamma+0.2 
$Cgam  =  cos($gamma*$pi/180) 

$b  =  $vol0/($a*$c*sqrt(l-$Calp*$Calp-$Cbet*$Cbet-$Cgam*$Cgam+2*$Calp*$Cbet*$Cgam)) 
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cellParameter  $a  $b  $c  $alpha  $beta  Sgamma 
writeFile  coordinate  filename  =  uc64-SbaNVT_def_$i.car 

#  Dynamics  Section:  NVT  equilibration 
minimize 
dynamics  \ 
time  =  10000.00  \ 
timestep  =  1.00  \ 

initial_temperature  =  298.00  \ 
ensemble  =  NVT  \ 
temperature  =  298.00  \ 
deviation  =  50000.00  \ 
execute  frequency  =  10.00  \ 

command  =  (print  table  file  =  uc64-SbaNVT_$i.tbl  +average  +batch_average  batch_size  =10  +cell 
+energy  +state  +stress  +strain)  \ 
execute  frequency  =  100.00  \ 

command  =  (print  archive  filename  =  uc64-SbaNVT_$i.arc  +coordinates)  \ 
execute  frequency  =  100.00  \ 

command  =  (writeFile  dynamics_restart  filename  =  uc64-SbaNVT_$i.xdyn) 
writeFile  coordinate  filename  =  uc64-SbaNVT_dyn_$i.car 


) 
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Appendix  C2 

%%  COMPILE  ALL  THE  DATA  CORRESPONDING  TO  THE  INDIVIDUAL  T RAJ. AC 
clc;  clear;  close  all; 

set (0, ’ Def aultFigureWindowStyle ' , 'docked ' ) ; 

natoms  =  68*64;  %  THIS  VALUE  MUST  BE  AJUSTED  FOR  EACH  SYSTEM 

endoftext  =  6537; 

nmbofrows  =  1000; 

nmbofcolms  =  51; 

avedata  =  zeros ( 70 ,  27 ) ; 

for  ii=l:75 

filename  =  sprintf ( ' . . /Sh64 J8Sac/uc64-SacNVT_%d . tbl ' , ii ) ; 
fid  =  f open ( filename ) ; 

tmpl  =  fscanf(fid,  %c ', endoftext ) ; 
tmpdata  =  zeros (nmbofrows , nmbofcolms ) ; 
data  =  zeros (nmbofrows, 26) ; 

%FROM  MATERIALS  STUDIO  [l.timestep  (even) .RunningAve  (odd) .BatchAve] 
%[l.timestep  2/3.TotalEnr  4/5.KinEnr  6/7.PotEnr  8/9. Temp  10/11. Press 
%12 /13 . Dens  14/15. Vol  16/17. a  18/19. b  20/21. c  22/23. alpha  24/25. beta 
%26/27. gamma  28/29. Sxx  30/31. Syy  32/33. Szz  34/35. Syz  36/37. Sxz 
%38/39.Sxy  40/41. exx  42/43. eyy  44/45. ezz  46/47. eyz  48/49. exz  50/51. exy] 
for  i=l : nmbofrows 

tmpdata (i,:)  =  f scanf ( f id,  ' %g ' ,  [  1  nmbofcolms]); 
end 

fid  =  fclose(fid); 

%Selecting  only  the  BATCH  data  (Odd  columns) 

% [l.timestep  2.TotalEnr  3.KinEnr  4 . PotEnr  5 . Temp  6. Press  7. Dens  8. Vol 
%9.a  10. b  11. c  12. alpha  13. beta  14. gamma  15. Sxx  16. Syy  17. Szz  18. Syz 
%19.Sxz  20.Sxy  21. exx  22. eyy  23. ezz  24. eyz  25. exz  26. exy  27.stn] 
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data (  :  ,  1 )  =  tmpdata (  :  ,  1 ) ; 
for  i=2:26 

data ( : , i)  =  tmpdata (:, 2*i-l ) ; 
end 

%Changing  the  sign  of  the  stress  values.  Accordint  to  AM 
data (:,  15:20)  =  (-1 . 0) *data ( : , 15 : 20) ; 

%Normalization  of  extensive  quantities  (Energies) 
data(:,2)  =  data ( : ,2) /natoms; 
data ( : , 3 )  =  data ( :  ,  3 ) /natoms ; 
data(:,4)  =  data (: ,  4 ) /natoms ; 

%Plots  for  the  energy  of  current  trajectory 

f igure ( ii ) /  plot (data ( :  ,  1 )  ,  data ( :  ,  2 )  ,  ' ko '  ,  1  Markers ize 1  ,  3 ) ; 
tmp  =  mean (data ( : , 2) ) /  axis([0  10e3  ( tmp+0 . 025^tmp)  ( tmp-0 . 025^tmp) ] ) ; 

xlabel('time  [fs]')/  ylabel ( sprintf (’ Energy  traj .  %d  ,ii)); 


% AVERAGES  FOR  EACH  DEFORMATION  STEP 
n  =  size (data, 1 ) ; 
if  (ii==l) 

avedata (ii, 1)  =  data (n,  1) /1000;  %Divide  by  1000  to  plot  in  [ps] 
else 

avedata (ii, 1 )  =  data (n,  1) /1000  +  avedata (ii-1 , 1 ) ; 
end 

avedata ( ii , 2 : 2 6 )  =  mean (data ( 101 : 1000 , 2 : 2 6) ) ; 
end 

closeall ; 
for  i=2 : 2  6 


figure (i) /  hold  on;  plot (avedata (:, 1 ), avedata (:, i ),’ ko MarkerSize ’ , 3 ) ; 


xlabel ( *  time  [ps]  ' ) /  ylabel ( sprintf ( ’ avedata ( %d)  , i ) ) ; 

switch  i 

case  2;  ylabel  ('E  [GPa] ’ ) ; 
case  8 

tmp  =  mean (data  (  : , 8 ) ) ;  axis([0  800  ( tmp-0 . 025*tmp)  ( tmp+0 . 025*tmp) 

ylabel ( ’ Volume ' ) ; 
case  15 

ylabel ( ’ Sxx  [GPa] ' ) ; 

save as (gcf ,  ’ Sxx-t ' ,  ' eps ’ ) ;  save as (gcf ,  ’ Sxx-t ' ,  ' png ’ ) ; 
case  16 

ylabel ( ’ Syy  [GPa]  ’  )  ; 

save as (gcf ,  ' Sy-t ' , ' eps 1 ) ;  save as (gcf, ! Syy-t ' , ' png ' ) / 
case  17 

ylabel ( ’ Szz  [GPa] ’ ) ; 

save as (gcf,  ’ Szz-t ’ ,  ' eps ' ) ;  save as (gcf,  ' Szz-t '  ,  ' png ’ ) ; 
case  18 

ylabel ( ' Syz  [GPa] ' ) ; 

save as (gcf, ’ Syz-t ' , ' eps ’ ) ;  save as (gcf, ' Syz-t ’ , ' png ' ) ; 
case  19 

ylabel ( ' Sxz  [GPa]  ?  )  ; 

save as (gcf,  ’ Sxz-t ' ,  ' eps ' ) ;  save as (gcf,  *  Sxz-t 1 ,  ' png ' ) ; 
case  20 

ylabel ( ’ Sxy  [GPa]  ' )  ; 

save as (gcf, ’ Sxy-t ’ , ' eps ’ ) ;  save as (gcf, ’ Sxy-t ' , ' png ' ) ; 

end 

end 


SHEAR  STRAIN  (ac) 
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avedata ( : , 27 )  =  (avedata (:, 13 ) -avedata ( 1 , 13 )) *pi/180 ; 

%%  PLOT  THE  DATA  OF  INTEREST 

%  FOR  THE  (beta, a)  SHEARING  SYSTEM:  tau(ac)  vs.  gamma (ac) 
figure  (21);  hold  on;  set (gca, 1 LineWidth ' , 2 . 0 ) ;  set (gca, ' FontSize  , 18 ) ; 
plot (avedata ( : , 27 ) , avedata ( :  ,  2 )  ,  ' ko ' , ’ Markers ize ’ , 4 , ’ MarkerFaceColor ’ , ’ k ’ ) ; 
axis  ( [0  0.3  -69.6  -69.3] ) ; 

%title (' Total  energy  vs.  Shear  strain,  \gamma_ (ac) FontSize ', 22 ) ; 

xlabel ( ’ \gamma_ (ac) ' , ' FontSize 1 , 22 , ' LineWidth ’,2.5) ; 

ylabel ( ' E  [ kcal*mol A ( -1 ) ] '  , ' FontSize ' , 22 , ' LineWidth ’,2.5) ; 

save as (gcf,  'uc64-Eac-stn',  ’ eps ’ ) ;  save as (gcf,  ’uc64-Eac-stn',  ’ png ' ) ; 

%  TRACTION  ON  a  PLANE 

alp  =  avedata ( 1 , 12 ) ;  bet  =  avedata ( 1 , 13 ) ;  gam  =  avedata ( 1 , 14 ) ; 
nor  =  [1  0  0] / 

drt  =  [cosd(bet)  cosd (alp) *sind (gam)  sqrt(l  -  (cosd (bet ) ) A2  - 
(cosd (alp) *sind (gam) ) A2 ) ] ; 
n  =  size (avedata (:, 1 ), 1 ) ; 
tra  =  zeros (n, 3);  % (x  y  z) 

%  from  avedata:  [15.Sxx  16.Syy  17.Szz  18.Syz  19.Sxz  20.Sxy] 
for  i=l:n 

tra(i,l)  =  avedata ( i , 15 ) *nor ( 1 )  +  avedata ( i , 2 0 ) *nor ( 2 )  + 

avedata (i, 19) *nor (3)  ; 

tra(i,2)  =  avedata ( i , 2 0 ) *nor ( 1 )  +  avedata ( i ,  1 6 ) *nor ( 2 )  + 

avedata (i, 18) *nor (3) ; 

tra(i,3)  =  avedata ( i , 1 9 ) *nor ( 1 )  +  avedata (i, 18) *nor (2)  + 
avedata (i, 17) *nor (3)  ; 


end 
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%  shear  (ac) :  component  of  the  traction  on  plane  a  along  c-dir:  (tra*c) 
shr  =  zeros (n, 1) / 
for  i=l:n 

shr(i)  =  tra (i, 1 ) *drt ( 1 )  +  tra (i, 2 ) *drt (2 )  +  tra (i, 3) *drt (3) ; 

end 

figure  (28);  hold  on;  set (gca ,  ' LineWidth ' ,  2 . 0 ) ;  set (gca,  ' FontSize ' , 18 ) ; 
plot (avedata ( : ,  27 ) ,  shr,  1 ko '  ,  '  Markers ize ' ,  4 ,  ' MarkerFaceColor ’ ,  ' k ’ ) ; 
axis  ([0  .3  min ( shr ) -0 . 25  max ( shr ) +0 . 25 ] ) ; 

%title (' Shear  stress  (ac)  vs.  Shear  strain  (ac)  ' ,  'FontSize'  ,  22) ; 

xlabel ( ’ \gamma_ (ac)  ,  ' FontSize '  ,  22  ,  LineWidth ' ,  2 . 5 )  ; 

ylabel ( ’ \tau_ (ac)  [GPa] ’ ,  ' FontSize ' ,  22 ,  ' LineWidth ' ,  2 . 5 ) ; 

save as (gcfA  ,uc64-Sac-stn,f  ’eps') ;  save as (gcfA  ,uc64-Sac-stn,f  1  png 1 ) ; 

figure  (29);  hold  on;  set (gca,  ’ LineWidth 2 . 0 ) ;  set (gca,  ' FontSize 1  ,  18 ) ; 
plot (avedata ( : , 27 )  ,  shr,  ko  *  j  ’ Markers ize ',4,  ' MarkerFaceColor ' ,  ' k ' ) ; 
plot (avedata ( : ,  27) ,  avedata ( : , 19) ,  'ko' , 'Marker Size ',4, 'MarkerFaceColor ' ,  ' y ' ) ; 
axis  ([0  .3  min ( shr ) -0 . 25  max ( shr ) +0 . 25 ] ) ; 

%title  (' Shear  stress  (ac)  vs.  Shear  strain  (ac)  ',' FontSize ', 22 )  ; 

xlabel ( ' \gamma_ (ac)  ' ,  ' FontSize '  ,  22 ,  ' LineWidth ',2.5)  ; 

ylabel (' \tau  [GPa]  ',' FontSize 1 , 22 ,' LineWidth 2 . 5 ) ; 

legend (gca, ' \tau_ (ac) ' , ' \tau_ (xz ) ' , ' Location ' , ' Northwest ' ) ; 

save as (gcf,  'uc64-SacSxz-stn',  'eps') ;  save as (gcf,  'uc64-SacSxz-stn',  ' png ' ) ; 


